You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
Can anyone tell me how to implement these timevarying state space equations in matlab
16 views (last 30 days)
Show older comments
Accepted Answer
Ced
on 22 Oct 2014
I would use one of the ode functions, e.g. ode45. You can then define an arbitrarily complex, nonlinear, time-varying, etc. system.
In your case, something like that (parameters are rubbish, but you have those)
function [Tout, Xout] = run_ode
% Some Parameters
Rm = 1000;
Ra = 1000;
% Time-varying parts:
C = @(t)(sqrt(t+1)^3);
dC = @(t)(3*sqrt(t+1)/2);
r = @(x)x*(-x<0); % Ramp function
p = @(x)[r(x(2)-x(1))/Rm ; r(x(1)-x(4))/Ra];
% Set up simulation
dx = @(t,x)ode_sys(t,x,C(t),dC(t),p(x));
Ts = 1e-12;
t_end = 1e-10;
t_sim = 0:Ts:t_end;
x0 = randn(5,1);
[ Tout, Xout ] = ode45(dx,t_sim,x0);
end
% Define system equations
function dy = ode_sys(t,x,C,dC,p)
% Some Parameters
Rs = 100;
Cr = 1e-6;
Cs = 2*Cr;
Ca = 0.1*Cr;
Ls = 6e-9;
Rc = 1e6;
b1 = 1/(Rs*Cr);
b2 = 1/(Rs*Cs);
% Defined time-varying system matrices elementwise
Ac = [ -dC/C 0 0 0 0 ;
0 -b1 b1 0 0 ;
0 b2 -b2 0 1/Cs;
0 0 0 0 -1/Ca;
0 0 -1/Ls 1/Ls -Rc/Ls ];
Pc = [ 1/C -1/C; -1/Cr 0; 0 0 ; 0 1/Ca ; 0 0 ];
dy = Ac*x + Pc*p;
end
Cheers
8 Comments
Bilal
on 22 Oct 2014
Thanks alot for all this, so here you have selected the C = @(t)(sqrt(t+1)^3); randomly ? secondly i need to extract x(1) x(2) x(3) x(4) x(5)
Bilal
on 22 Oct 2014
this is how i have calculated the time varying C(t) and dC(t), can you please tell me how to write C(t) which is the reciprocal of E(t) in ode45 syntax. C(t)=1/E(t)
clc;
close all;
clear all;
%%
Emax = 2;
Emin = 0.06;
HR = 60;
ts = 0.0025;
t=[0:ts:1];
tc = 60/HR;
Tmax = 0.2 + 0.15 *tc;
tn=t/Tmax;
En= 1.55*(((tn./0.7).^1.9)./(1+(tn./0.7).^1.9)).*((1./(1+(tn./1.17).^21.9)));
E= ((Emax - Emin) * En)+Emin;
plot(t,E);
grid on
%% for dC/dt
C=[1./E];
dcdt=diff(C)./ts;
figure
C=C(1:length(C)-1);
CC=dcdt./C;
plot(t(1:(length(t)-1)),-CC);
title('\Delta{C}/\Delta{t}')
grid on
Ced
on 22 Oct 2014
Yes, I chose C(t) and dC(t) randomly.
One important thing is not quite clear to me. Do your functions C,dC, and E have to be discrete? or can you define them as continuous functions of time t?
Ced
on 22 Oct 2014
Edited: Ced
on 22 Oct 2014
Assuming you can define them as continuous functions: You don't need any particular syntax for C(t), dC(t) with ode45. The only function which needs a particular syntax is the function "dx" which I passed to ode45 itself. "dx" needs to have the variables (t,x). For the rest, you just need to define "standard" matlab functions.
Basically, there are two ways to define functions:
a) anonymous functions: "inline" functions as above with the @ sign
b) functions defined outside of the rest of your code. The syntax there would be "function output = function_name(input)". You can have them as separate .m files, or locally/nested in your main script file. All details are in the link, but this is not so important for the problem at hand.
You can find a good introduction with examples etc here: http://www.mathworks.com/help/matlab/functions.html
In your example, I would say that anonymous functions are enough, but I don't know the full complexity of your program. If there are a lot of interdependent and/or complex functions, I would define them as "local functions". See the link above for examples.
If you want to define E(t) as an anonymous function:
En = @(tn)1.55*(((tn./0.7).^1.9)./(1+(tn./0.7).^1.9)).*((1./(1+(tn./1.17).^21.9))); % note that this is a function, and not a vector of numbers
E = @(t) ((Emax - Emin) * En(t/Tmax))+Emin; % again, E(t) as a function
Let me know if something is unclear.
Bilal
on 22 Oct 2014
so do i need to give the range of tn ??
and what if i want to take the derivative of E here as i need C and dC which is C=1/E so from above equation E do i need to manually calculate the derivative can't i use the diff function any more as its a function
Ced
on 22 Oct 2014
no, you don't need to give the range of tn for ode45. You just define E the same way you do for C, dC and the rest. When you then call e.g. E(3), then the function E(t) will be evaluated at t = 3;
The integration range is then passed to ode45. In my first answer, I had used
Ts = 1e-12;
t_end = 1e-10;
t_sim = 0:Ts:t_end;
x0 = randn(5,1);
[ Tout, Xout ] = ode45(dx,t_sim,x0);
t_sim is only passed to ode45, and does not appear anywhere else. Instead of defining t_sim as a time vector with all time steps, I could also only pass [ t_start t_end ], and then ode45 output it's own time steps. The reason I prefer passing t_sim as "full vector" is that it makes things easier to plot later, but that is a personal preference.
For the derivative: Yes, calculate the derivative by hand, and then hardcode the function. There are ways to compute derivatives symbolically using the symbolic math toolbox, but that would just make things unnecessarily complicated here. I think this derivative is still ok. If you want to check that your derivative is correct, you can always use some calculator or mathematica to check your answer. Also, note that "diff" does not compute derivatives as such, but only takes the difference between the elements in your vector, e.g. diff([ 1 2 3]) would return [ 1 1 ].
More Answers (6)
Bilal
on 22 Oct 2014
yes you are right but the thing is E(t) is the re scaled version of En(tn) so it would definitely effect E(t) if i don't define tn anywhere if you run my code, you will see i m getting the write results as given in this paper. check the attachment
And yes if i accept the answer can we still continue communicating here. ?
26 Comments
Ced
on 22 Oct 2014
Edited: Ced
on 22 Oct 2014
Yes, I'm pretty sure we can keep adding to the answer once you accept it. Worst case, email me afterwards.
Thanks for the paper. But you do add Emin, Emax, and tc as known constants, correct? So tn is nothing else than t/Tmax, meaning that instead of writing E(t) = (Emax-Emin)*En(tn)+Emin, write it as a function of t, i.e.
E = @(t) ((Emax - Emin) * En(t/Tmax))+Emin;
Of course, if you really want, you could define
tn =@(t)t/Tmax
as well an use it like this, but I see no reason to do this. The only reason they introduce tn is that it makes the formula of En(tn) more compact for the paper.
*EDIT* Just realized this may not have come across clearly: We are only talking about ode45 here. If you want to plot the functions, you of course have to tell 'plot' at which points you want your functions evaluated. But let's maybe stick to solving the ode first, and then we can tackle the plotting part.
Bilal
on 22 Oct 2014
Edited: Bilal
on 22 Oct 2014
i don't have any problem in writing t/Tmax or tn both are same all i wanted to know if i write express like this E = @(t) ((Emax - Emin) * En(t/Tmax))+Emin; in ode45 it means i need to have the values of all the contants except (t) as the function is dependent on t. Am i right, yes everything is known Emin Emax and tc.
Bilal
on 22 Oct 2014
I have modified the code but now i am getting a very strange error in ode45. please have a look
Ced
on 22 Oct 2014
Yes, the error message is not so intuitive. The problem is that in your definitions of En(t), C(t), etc., which are functions of "t", you are using "tn". You did define "tn" as a function, which is fine, but this means that you have to pass an argument when you use it. In short: just replace "tn" by "tn(t)" in your definitions of En, E, C, etc., and you're good to go.
I'll have to pass for today, but I'll be back tomorrow.
Cheers
Bilal
on 22 Oct 2014
Thanks alot for your effort and time it worked but u know what now Xout matrix is undefined having (NaN) in it
Ced
on 23 Oct 2014
Yes, that's unfortunate :D. My best guess at the moment: Your initial conditions are still set as "randn" as far as I can tell. 'randn' takes random values from the normal distribution. If nothing is known, this is usually as good a guess as any other, but since you have a biological system, there may be some limitations to what makes sense.
Things you could do:
- Check all equations (in particular, signs)
- Check if the initial conditions are given somewhere in the paper.
- check if there are any upper/lower bounds on the values of x (e.g. a negative pressure makes little sense to me), make sure the initial conditions fulfil them
- There is most likely a division by zero somewhere. Try to find out where. To do this, you can e.g. write 'keyboard' in your code. The run will then pause there, and you can evaluate the code one line at a time by pressing f10. There are lots of other debugging tools for matlab, you can check the documentation for a list.
- your time scale seems to be similar to the ones in the paper, so that's good, and since the curve is very smooth and does not have high frequency noise or so, I don't think it's the solver tolerance. You could however try another solver (look at 'doc ode45' for what other versions are available), but since your solution is already NaN after 1 timestep, I don't think it's the solver.
Bilal
on 23 Oct 2014
yes there was something wrong in my dC equation (the derivative of c equation) i choose dc = C just to check the NaN problem it worked but now all the values are same means nothing is updating, and there is no info in the paper regarding initial condition but if you see the plots all the values are positive so i made randn to rand but the values are small is that okay as its just the initial condition. one thing i wanna ask here is there any way to plot E(t) here, so that i can verify that i m getting the same plot as i was getting before and then i can plot C and dC if i m getting the same plots then I would shift to next step. Hope you can understand what i mean to say
Ced
on 23 Oct 2014
Glad it works! Not sure I understand, no. What do you mean by "plot here"? Since E, C, and dC are only functions of time, they are not dependent on your solution from the ode, so I'm not sure what exactly you want to check.
But you can plot them whenever you like using "plot", e.g.
plot(t_sim,E(t_sim))
Bilal
on 23 Oct 2014
yup i did it and i m getting the same plot
function [Tout, Xout] = run_ode
% Some Parameters
Rm = 0.0050;
Ra = 0.0010;
Emax = 2;
Emin = 0.06;
HR = 60;
tc = 60/HR;
Tmax = 0.2 + 0.15 *tc;
ts=0.0025
t=[0:ts:1];
% Time-varying parts:
tn = @(t) t./Tmax;
En = @(t) 1.55*(((tn(t)./0.7).^1.9)./(1+(tn(t)./0.7).^1.9)).*((1./(1+(tn(t)./1.17).^21.9)));
E = @(t) ((Emax - Emin) * En(tn(t)))+Emin;
plot(t,E(t))
Bilal
on 23 Oct 2014
but i don't know how to take the derivative of this complex function E, i tried online solvers for derivative i am getting very wired expression from there.. will solve it then will get back to you.
Ced
on 24 Oct 2014
Edited: Ced
on 24 Oct 2014
Sounds good. Yes, solutions can look strange sometimes because they apply weird "simplifications". You could also simply derive it by hand, the first derivative is still fairly simple, and then compare the plot to the one from your online solution. Let me know if you need help with that.
Bilal
on 25 Oct 2014
can you please check the problem in the code, getting some sort of error in ode_sys. Please have a look at the attached file.
Ced
on 26 Oct 2014
Edited: Ced
on 26 Oct 2014
ok, had a look. A new version is in the attachments, but a couple of (hopefully) useful comments:
The error you are (were) getting is:
Error using vertcat Dimensions of matrices being concatenated are not consistent.
Error in tmp>ode_sys (line 60) Ac = [ -dC/C 0 0 0 0 ;
Error in tmp>@(t,x)ode_sys(t,x,C(t),dC(t),p(x))
[ ... ]
Error in tmp (line 45) [ Tout, Xout ] = ode45(dx,t_sim,x0);
Generally speaking, the first line tells you what kind of error it is, and then goes on telling you where this error was produced. In this case, the error is that "dimensions of matrices being concatenated are not consistent", and it happens in (line 60).
This is quite a common error which means that you are trying to concatenate two matrices/vectors (meaning "insert" them into the same array), but that there dimensions don't match. This usually happens when you messed up indices somewhere in a for loop or so. Unfortunately, in your case, there's a bigger problem. What you are doing with "diff" is that you are approximating the derivative by finite differences. But ode45 evaluates your function at a single timestep, and with only one value of C given, computing a finite difference is impossible.
I somewhat fixed your code with a hack which I think should work, but I would really encourage you to use the "real" derivative. I'll see if I find some time tomorrow and set it to you.
A couple of useful side notes:
1) The last index can be reach through "end", e.g.
x(1:end-1) % x without last element
which is a bit easier than your t(1:length(t)-1) construct :).
2) I wouldn't implement a function C referring to a previous version of C. I'm actually surprised this even worked. I simply omitted dC and plugged C directly into CC. If for some reason, you ever had to do that, I would use two different names, e.g. C0 and C1
Cheers
Bilal
on 27 Oct 2014
Oh thanks alot it worked, now coming back to the state space, if you see in the paper, there are five state variables, x1 to x5. they r in dy in our code if im not getting it wrong, but i m getting only single value of x1 x2 x3 x4 x5. I need a vector so that i can plot these values.
Ced
on 27 Oct 2014
Hi
Sorry it took me so long. I wanted to check our equations, but am running into big problems, nothing seems to match. Probably some stupid typo somewhere. Anyway... to your question: Tout is a vector of time, Xout is a matrix containing your states. So x1, x2, ... are in Xout. dy in our code is x_dot.
if your run the function from the matlab command line (you need to comment out "clear all" first in the function, otherwise it won't work), each xi will be in one column, i.e.:
[ t_out, x_out ] = run_ode;
x1 = x_out(:,1);
x2 = x_out(:,2);
etc.
Bilal
on 27 Oct 2014
I know but x1 x2 ... n x5 never gets updated their values remain same all the time.
Bilal
on 27 Oct 2014
although if you see C(t) and dC(t) is changing but the state variable never changes.
Ced
on 28 Oct 2014
Ah, yes. Two things:
1. When computing the derivative with finite differences, you need to take the same "step size" on the top and bottom, e.g. for f([t t+dt1])/dt2, both dt1 and dt2 have to be the same (my mistake, sorry).
Replace the following lines in your code, this should solve the problem:
CC = @(t,ts) diff(C(t))./(C(t(1:end-1))*ts); % dC./C
[...]
plot(t(1:end-1),-CC(t,ts));
[...]
Ts = 0.0025;
dx = @(t,x)ode_sys(t,x,C(t),CC([t t+Ts],Ts),p(x)); % Reverse order!
2. Then, concerning the non-changing data: Have a look at your time step Ts and the end time. You are only simulating for 1e-10 seconds! During this time interval, nothing changes. You could e.g. set Ts = 0.0025 and t_end = 1 for one cycle, but that depends on what you want.
Bilal
on 28 Oct 2014
r = @(x)x*(-x<0); % Ramp function p = @(x)[r(x(2)-x(1))/Rm ; r(x(1)-x(4))/Ra];
how this statement will get the values of x(1) x(2) and x(3) it is written above the ode_sys call ... i m not getting it.
Bilal
on 28 Oct 2014
I have tried to solve the problem, now its the update version of code pls run it and check that it only calculates states for 1 cycle. although now i have 4 or 5 cycles of time varying signal. so the plot in Xout (:,1) ... Xout(:,5) should repeact periodically like E(t) and C(t) ...
And attached herewith are the plots from paper that i have already sent you.
Ced
on 28 Oct 2014
I'm sorry, I don't understand your question(s). Why don't you just plot it? That way, you'll see how many cycles appear? So... just to e.g.
plot(Tout,Xout(3,:))
And if it's not enough/too many cycles, adapt your end time accordingly (or simply only plot a particular section of the results).
Bilal
on 29 Oct 2014
4 Comments
Ced
on 29 Oct 2014
Yes, of course. Since you know Xout, you can also compute p, r, and everything else. If you want to have them as "standard outputs" like Xout, you could do something like :
Nt = length(Tout);
Pout = zeros(Nt,2);
Rout = r(Xout);
for i = 1:Nt
Pout(i,:) = p(Xout(i,:));
end
after Tout and Xout have been computed and add them to the output, i.e.
function [Tout, Xout, Pout, Rout] = run_odebq
IMPORTANT: You will need to change your ramp function to an element-wise evaluation, otherwise you get some errors. Meaning it should look like this (note the '.'):
r = @(x)x.*(-x<0); % Ramp function
Once you have Pout, Rout, you can plot them, do some analysis or whatever (you can of course also plot them directly in the function without passing them as outputs)
Bilal
on 29 Oct 2014
okay, i got it thanks alot.
Did you find anything why i am getting only one cycle in Xout and then output gradually goes to zero
Ced
on 30 Oct 2014
No, I'm afraid I'm a bit caught up in work myself at the moment. I'll have a look at it if I find the time. I haven't read the paper in detail, but maybe they have something like mod(t,HR) or so? Maybe you can try to find out what function is not periodic, and if not, why.
Bilal
on 31 Oct 2014
This is what exactly i have done in the generation of C(t) i have use mod(t,HR) but do i need to implement it for the state space as well ? okay please get back to me when u have time. thanks
Ced
on 31 Oct 2014
I mean, that would certainly work. simply replace every instance of "t" by "mod(t,60/HR)" when calling ode_sys, so
dx = @(t,x)ode_sys(mod(t,60/HR),x,C(mod(t,60/HR)), ...
The question is if this makes sense. That's something you could ask your professor or discuss with some phd student. Skimming the paper, I cannot find any reference that they artificially periodicized (that's a word now ^^) their solutions. But then again, little details are often omitted (or forgotten) when describing simulations in papers.
12 Comments
Bilal
on 31 Oct 2014
Thanks alot by doing this the periodicity issue has been resolved, I am getting the shape of all the states exactly like the results on paper (Except X3).
Bilal
on 31 Oct 2014
r = @(x)x*(-x<0); % Ramp function
Can you please tell me how this is interpreted.
we need output equal x as long as x>0 so how the above statement works
Bilal
on 31 Oct 2014
Secondly if you see my plots they change on every cycle, whereas in paper there is no change, the plot is periodic.
Ced
on 2 Nov 2014
Great!
What you want the ramp function r(x) to do is the following: If the value x is greater or equal than zero, set r(x) = x, else set r(x) = 0. The
(-x<0)
part returns a logical value, which is 1 if (0 <= x) and zero else (if used in an algebraic expression, the logical is automatically interpreted as a double). I wrote it as -x<0 to have an inequality instead of a <=.
So, x*(-x<0) is x*1 if 0 <= x and 0 else, which is just what the ramp function does.
Now, to the changing period: This is a little strange. To be fair, the first period in the paper is also slightly different than the others. Check if the mod(...) arguments are the same everywhere in your code. Unfortunately, it could be that the equations are very sensitive to changes in the derivative, and since our dC is approximated, that could explain the differences. You could test this by changing the magnitude of the step in the finite difference computation (the second argument of dC). If that does not change anything, let me know, and I'll think a bit harder :)
Ced
on 3 Nov 2014
hm... looks fishy to me, but I can't really put my finger on the reason. One thing I would certainly change is that you should not use randn for your initial conditions, since your pressure should not be negative.
And then, although I don't think that the derivative is really the problem, there is a strong dependency on it, as the solution changes when I change the discretization. That should not happen. If you could find the algebraic derivative, that would certainly help a lot. If you can't figure out the correct derivative by hand, you could try using the symbolic math toolbox, e.g.
syms t_smb real % define symbolic variables
f = t_smb.^2 + 3; % define symbolic expression (would be your function)
df = diff(f,t_smb); % compute derivative wrt t_smb
df_fun = matlabFunction(df); % transform into matlab function
Note that here, "diff" is overloaded and computes the algebraic derivative and not the finite difference, since your input is symbolic.
Also, I actually find it a bit weird to have mod(t,60/HR). It would make more sense to have mod(t,Tmax) since t is "normalized" by Tmax. However, this diverges. Do you have anyone you can ask (supervisor or so) about this periodicity and normalization who is working in that field? While I know how ode works, I don't really know anything about heart simulations.
Bilal
on 4 Nov 2014
yes i will ask my supervisor. getting too wired expression through this method pls check did i write the expression correctly ?
if true
% code
end
Bilal
on 7 Nov 2014
I am unable to plot derivative of C that you told me to generate using syms. its getting mpower error
See Also
Categories
Find more on Startup and Shutdown in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!An Error Occurred
Unable to complete the action because of changes made to the page. Reload the page to see its updated state.
Select a Web Site
Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .
You can also select a web site from the following list
How to Get Best Site Performance
Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.
Americas
- América Latina (Español)
- Canada (English)
- United States (English)
Europe
- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)
- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)
Asia Pacific
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)