You are now following this question
- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.
I need help solving a system of differential equations. The equations are given below, in matrix form. The problem that I'm having is regarding the fact that I have time dependant elements in the matrices.
2 views (last 30 days)
Show older comments
42 Comments
Walter Roberson
on 4 Aug 2020
That should not be a problem if and are known. Just write the matrix equation as-is, replacing and with variables, and then insert a line or two to assign the variables based upon the current time input parameter and the known functions.
Jelena Kresoja
on 4 Aug 2020
Yes, thank you! This is what I initially wanted to do but I'm not sure how to assign the variable based on time.
James Tursa
on 4 Aug 2020
What do your C(t) and Cdot(t) functions look like? An example would be
function dx = myderivative(t,x)
C = t^2 + 1;
Cdot = 2*t;
dx = all of your matrix stuff above
end
Jelena Kresoja
on 5 Aug 2020
So, C=1/E, and Cdot I also have to compute. In the eqution below tn=t/Tmax where Tmax is know.
Walter Roberson
on 5 Aug 2020
I do not understand the notation with regards to implying that E is a vector whose minima is to be taken, but implying that instead E is a family of functions ??
Walter Roberson
on 5 Aug 2020
If your E is a vector, and is min(E) and is max(E), then what does mean? With E being a vector, I could see being notation for the element of E . Does mean the element of E, multiplied by the element of t ?
Sam Chak
on 6 Aug 2020
If is a scaled quantity of time, and is continuously differentiable, then is a continuous function. I don't know the value of , but the pattern of should look similar to this:
t = linspace(0, 2, 1001);
Cdot = (0.224843*t.^24.7 + 0.104268*t.^22.8 - 0.308422*t.^0.9)./(0.0006283*t.^23.8 + 0.000319046*t.^21.9 + t.^1.9 + 0.00993399).^2;
plot(t, Cdot)
Jelena Kresoja
on 6 Aug 2020
Yes, Sam, thats how its supposed to look. Tmax=0.2+0.15*tc, where tc=60/heart_rate. Im using 60 beats per minute for heart rate. And when i have all of that i just dont know how to put that in the matrix so i use just the value at a given moment?
Sam Chak
on 6 Aug 2020
You can follow the advices from Walter Roberson and James Tursa, or refer to the examples in this link:
Here is a simple example to simulate the system when there are time-varying components:
The codes
clear all; clc
tspan = 0:0.01:2;
y0 = [1];
[t, y] = ode45(@(t, y) [(-((0.224843*t^24.7 + 0.104268*t^22.8 - 0.308422*t^0.9)/(0.0006283*t^23.8 + 0.000319046*t^21.9 + t^1.9 + 0.00993399)^2)/(1/((2 - 0.06)*(1.55*(((t/0.7)^1.9)/(1 + (t/0.7)^1.9))*(1/(1 + (t/1.17)^21.9))) + 0.06)))*y(1)], tspan, y0);
% Cdot/C
CC = ((0.224843*t.^24.7 + 0.104268*t.^22.8 - 0.308422*t.^0.9)./(0.0006283*t.^23.8 + 0.000319046*t.^21.9 + t.^1.9 + 0.00993399).^2)./(1./((2 - 0.06)*(1.55*(((t/0.7).^1.9)./(1 + (t/0.7).^1.9)).*(1./(1 + (t/1.17).^21.9))) + 0.06));
subplot(2,1,1)
plot(t, y)
subplot(2,1,2)
plot(t, CC)
Results
Explanation
Because this is a 1st-order system, when , will diverge, and when , will converge to a steady-state value.
Jelena Kresoja
on 8 Aug 2020
Can I solve my initial system of equations in matrix form using dsolve? If yes could I please just get an example of how do use it in a situation when I have the above mentioned time dependant elements in the matrix? This is what I dont know how to do right now.
Walter Roberson
on 8 Aug 2020
I can't tell at the moment whether dsolve would work; I would have to type all those equations in to try it.
dsolve() does not inherently have problems with time-varying terms. However, when I look at your it looks messy enough to me that I would not expect dsolve() to be able to find closed form solutions.
The symbolic toolbox can still be useful, though, in that it provides a way to code the equations in more understandable form. Then you follow the workflow show in the first example of odefunction() to convert into something that can be processed by ode45()
Jelena Kresoja
on 9 Aug 2020
My problem with ode45 is that I dont undrestand, for my case, how to define the ode, because for example x1 is dependant on x2 and x4 which also need to be integrted. And I dont quite understand help thats written for this. I would try what James Tursa wrote in one of the previous comments:
function dx = myderivative(t,x)
C = t^2 + 1;
Cdot = 2*t;
dx = all of your matrix stuff above
end
i just dot know in which form should I write all of the elements of the matrix?
Walter Roberson
on 9 Aug 2020
I recommend using the symbolic toolbox.
syms t
syms x1(t) x2(t) x3(t) x4(t) x5(5)
x = [x1; x2; x3; x4; x5];
C(t) = some expression
Cdot = diff(C,t);
M1 = [-Cdot./C, 0, 0, 0, 0, 0;
0, -1./(Rs.*Cr), 1./(Rs.*Cr), 0, 0, 0;
and so on]
M2 = [1./C, -1./C;
-1./Cr, 0;
and so on]
M3 = [r.*(x(2)-x(1))./RM;
r.*(x(1)-x(4))./Ra];
dx = diff(x);
eqn = dx == M1 * x + M2 * M3;
Now follow the work flow at the first example of odeFunction() in order to turn eqn into something that can be processed numerically by an ode*() routine.
Jelena Kresoja
on 10 Aug 2020
Im getting an error message:
Error using symengine
Dimensions do not match.
The problem is with matrices M2 ad M3, specifically (x(2)-x(1)), without this it works fine.
Walter Roberson
on 10 Aug 2020
Edited: Walter Roberson
on 10 Aug 2020
What is size(RM) and size(Ra) ?
Is the mismatch error in creating M3, or is it in M2 * M3 in building eqn ? If it is M2 * M3 then what is size(M2) and size(M3) ?
Jelena Kresoja
on 10 Aug 2020
Ra and Rm are just constants. The error is for M2*M3. The size of 1x1 and for M3 its 10x1.
Walter Roberson
on 10 Aug 2020
M2 should not be 1 x 1: it should be 5 x 2. And M3 should be 2 x 1. Unless your r is 5 x 1 ??
Note that I read your equations as r being a constant being multiplied by x(2)-x(1), rather than r being a function being called with parameter x(2)-x(1)
Jelena Kresoja
on 11 Aug 2020
Oh the matrices themselves are these dimention. M2 is 5x2 and M3 is 2x1. I also tried to put in the function r but it doesnt work. Even if I just write x(2)-(1) I get an error message.
Jelena Kresoja
on 11 Aug 2020
Rs = 1; %systemic resistance
Rm = 0.005; %mitral valve resistance
Ra = 0.001; %aortic valve resistance
Rc = 0.0398; %chracteristic resistance
Cr = 4.4; %left atrial compliance
Cs = 1.33; %systemic compliance
Ca = 0.08; %aortic compliance
Ls = 0.0005; %inertance in aorta
hr = 60; %hear rate
Emax = 2; %max elastance
Emin = 0.06; %min elastance
tc = 60/hr;
t = 0:0.01:tc;
Tmax = 0.2+0.15*tc;
syms t
syms x1(t) x2(t) x3(t) x4(t) x5(t)
x = [x1; x2; x3; x4; x5];
C(t) = 1./((Emax-Emin)*(1.55.*(((((t./Tmax)./0.7).^1.9)./(1+(((t./Tmax)./0.7).^1.9))).*(1./(1+(((t./Tmax)./1.17).^21.9)))))+Emin);
Cdot = diff(C,t);
M1 = [-Cdot./C, 0, 0, 0, 0;
0, -1./(Rs.*Cr), 1./(Rs.*Cr), 0, 0;
0, 1./(Rs.*Cs), -1./(Rs.*Cr), 0, 1./Cs;
0, 0, 0, 0, -1./Ca;
0, 0, -1./Ls, 1./Ls, -Rc./Ls];
M2 = [1./C, -1./C;
-1./Cr, 0;
0, 0;
0, 1./Ca;
0, 0];
M3 = [((x(2)-x(1)).*(x(2)-x(1)>=0))./Rm;
((x(1)-x(4)).*(x(1)-x(4)>=0))./Ra];
dx = diff(x);
eqn = dx == M1*x + M2*M3;
Walter Roberson
on 11 Aug 2020
please format as code. I cannot copy and paste that from my mobile phone.
Jelena Kresoja
on 11 Aug 2020
Rs = 1; %systemic resistance
Rm = 0.005; %mitral valve resistance
Ra = 0.001; %aortic valve resistance
Rc = 0.0398; %chracteristic resistance
Cr = 4.4; %left atrial compliance
Cs = 1.33; %systemic compliance
Ca = 0.08; %aortic compliance
Ls = 0.0005; %inertance in aorta
hr = 60; %hear rate
Emax = 2; %max elastance
Emin = 0.06; %min elastance
tc = 60/hr;
t = 0:0.01:tc;
Tmax = 0.2+0.15*tc;
syms t
syms x1(t) x2(t) x3(t) x4(t) x5(t)
x = [x1; x2; x3; x4; x5];
C(t) = 1./((Emax-Emin)*(1.55.*(((((t./Tmax)./0.7).^1.9)./(1+(((t./Tmax)./0.7).^1.9))).*(1./(1+(((t./Tmax)./1.17).^21.9)))))+Emin);
Cdot = diff(C,t);
M1 = [-Cdot./C, 0, 0, 0, 0;
0, -1./(Rs.*Cr), 1./(Rs.*Cr), 0, 0;
0, 1./(Rs.*Cs), -1./(Rs.*Cr), 0, 1./Cs;
0, 0, 0, 0, -1./Ca;
0, 0, -1./Ls, 1./Ls, -Rc./Ls];
M2 = [1./C, -1./C;
-1./Cr, 0;
0, 0;
0, 1./Ca;
0, 0];
M3 = [(x(2)-x(1))./Rm; (x(1)-x(4))./Ra];
dx = diff(x);
eqn = dx == M1*x + M2*M3;
Walter Roberson
on 11 Aug 2020
M2 = [1./C(t) , -1./C(t);
-1./Cr, 0;
0, 0;
0, 1./Ca;
0, 0];
M3 = [(x2(t)-x1(t))./Rm; (x1(t)-x4(t))./Ra];
dx = diff(x);
eqn = dx == M1(t)*x(t) + M2*M3;
Walter Roberson
on 11 Aug 2020
eqn will look a bit messy. That is due to matlab converting those floating point powers into rationals and expanding out numerators. There is a way to prevent that.
Walter Roberson
on 11 Aug 2020
No error message for me.
Rs = 1; %systemic resistance
Rm = 0.005; %mitral valve resistance
Ra = 0.001; %aortic valve resistance
Rc = 0.0398; %chracteristic resistance
Cr = 4.4; %left atrial compliance
Cs = 1.33; %systemic compliance
Ca = 0.08; %aortic compliance
Ls = 0.0005; %inertance in aorta
hr = 60; %hear rate
Emax = 2; %max elastance
Emin = 0.06; %min elastance
tc = 60/hr;
t = 0:0.01:tc;
Tmax = 0.2+0.15*tc;
syms t
syms x1(t) x2(t) x3(t) x4(t) x5(t)
x = [x1; x2; x3; x4; x5];
C(t) = 1./((Emax-Emin)*(1.55.*(((((t./Tmax)./0.7).^1.9)./(1+(((t./Tmax)./0.7).^1.9))).*(1./(1+(((t./Tmax)./1.17).^21.9)))))+Emin);
Cdot = diff(C,t);
M1 = [-Cdot./C, 0, 0, 0, 0;
0, -1./(Rs.*Cr), 1./(Rs.*Cr), 0, 0;
0, 1./(Rs.*Cs), -1./(Rs.*Cr), 0, 1./Cs;
0, 0, 0, 0, -1./Ca;
0, 0, -1./Ls, 1./Ls, -Rc./Ls];
M2 = [1./C(t) , -1./C(t);
-1./Cr, 0;
0, 0;
0, 1./Ca;
0, 0];
M3 = [(x2(t)-x1(t))./Rm; (x1(t)-x4(t))./Ra];
dx = diff(x);
eqn = dx(t) == M1(t)*x(t) + M2*M3;
%now follow odeFunction first example
[eqs,vars] = reduceDifferentialOrder(eqn,x(t));
[M,F] = massMatrixForm(eqs,vars);
f = M\F;
odefun = odeFunction(f,vars);
%hen
initConditions = [vector of 5 values]; %[x1(0), x2(0), x3(0), x4(0), x5(0)];
tspan = [0 10]; %or as appropriate
%and run
[t, x] = ode45(odefun, tspan, initConditions);
Jelena Kresoja
on 11 Aug 2020
This works for me too now....
I just still have to change the elements of matrix M3, and if I try and do that it doesnt work again. The function that is actually r is 0 when x2(t)-x1(t) is less than 0 and keeps the value x2(t)-x1(t) when this is larger than 0.
Walter Roberson
on 11 Aug 2020
You must use event functions for that situation, as it is not differentiateable at the boundary.
Walter Roberson
on 12 Aug 2020
In MATLAB in numeric mode,
r = @(expression) expression.*(expression > 0);
unless expression can be infinite. If the expression can be infinite then more work has to be done.
The version with sign() fails for negative infinity.
In MATLAB in numeric mode, you can also use
r = @(expression) max(0, expression)
which does not have the problem with -infinity. However, I tend to avoid using min() and max() for this kind of work as it is not uncommon for me to want to be able to pass in symbolic expressions, and min() and max() do not play nicely with symbolic expressions.
All of these versions, using sign() or .*(logical) or max() or min(), cannot be differentiated at 0. The same is also true if you code using piecewise(): you cannot differentiate at 0...
... And that is important for the ode*() functions. The ode*() functions all require that the first derivative of the coded functions are continuous, as the mathematics that they are designed against for optimal evaluation points require first derivatives to be continuous.
Therefore when you want to switch between 0 and not-zero, you must use an event function to stop integration, and then you restart integration.
Jelena Kresoja
on 12 Aug 2020
So i wrote this is as:
r = @(ksi) ksi.*(ksi>0);
options = odeset('Events',r);
And then I just called the r function where needed (as r(x2(t)-x1(t))), but then I get an error. I probably missed some steps?
Jelena Kresoja
on 13 Aug 2020
Okay, so I made an event function:
function [val, isterminal, dir] = eventsFun(t,x)
val = x(2)-x(1);
isterminal = 0;
dir = 0;
end
and then i just added it to odeset.
Is this correct now or do I need something more? Cause my results arent any different.
Walter Roberson
on 13 Aug 2020
If the results are not any different then possibly the event was never triggered.
You can use the 5-output version of ode45() to get information about the event functions that were triggered.
Jelena Kresoja
on 14 Aug 2020
So the events were triggered but nothing happened I guess?
I think im missing something cause the event function is detecting when x2-x1=0 but I need the value to be set to 0 if x2-x1 is also less than 0. If its greater than zero it needs to stay x2-x1.
Walter Roberson
on 14 Aug 2020
You have to loop on the ode45() call, like
mint = 0; maxt = 60;
x0 = appropriate vector
all_t = {}:
all_x = {};
while true
[t, x] = ode45(YourFunction, [mint maxt], x0, opts);
all_t{end+1} = t;
all_x{end+1} = x;
mint = t(end);
if mint == maxt; break; end %we are done
x0 = x(end,:);
end
It is okay if YourFunction has if statements in it, as long as they always evaluate to the same thing for any one call to ode45()
Jelena Kresoja
on 14 Aug 2020
Im so sorry i literally dont get whats written here, or where im supposed to put this.....
Answers (2)
hosein Javan
on 10 Aug 2020
if you are solving a circuit with time-variant capacitors then your state equations are no longer in "Xdot=A*x+B*U", and they are expressed in general form "Xdot=f(X,U)". you have use ode solvers.
Walter Roberson
on 14 Aug 2020
Rs = 1; %systemic resistance
Rm = 0.005; %mitral valve resistance
Ra = 0.001; %aortic valve resistance
Rc = 0.0398; %chracteristic resistance
Cr = 4.4; %left atrial compliance
Cs = 1.33; %systemic compliance
Ca = 0.08; %aortic compliance
Ls = 0.0005; %inertance in aorta
hr = 60; %hear rate
Emax = 2; %max elastance
Emin = 0.06; %min elastance
tc = 60/hr;
t = 0:0.01:tc;
Tmax = 0.2+0.15*tc;
syms t
syms x1(t) x2(t) x3(t) x4(t) x5(t)
x = [x1; x2; x3; x4; x5];
r = @(v) heaviside(v) * v
C(t) = 1./((Emax-Emin)*(1.55.*(((((t./Tmax)./0.7).^1.9)./(1+(((t./Tmax)./0.7).^1.9))).*(1./(1+(((t./Tmax)./1.17).^21.9)))))+Emin);
Cdot = diff(C,t);
M1 = [-Cdot./C, 0, 0, 0, 0;
0, -1./(Rs.*Cr), 1./(Rs.*Cr), 0, 0;
0, 1./(Rs.*Cs), -1./(Rs.*Cr), 0, 1./Cs;
0, 0, 0, 0, -1./Ca;
0, 0, -1./Ls, 1./Ls, -Rc./Ls];
M2 = [1./C(t) , -1./C(t);
-1./Cr, 0;
0, 0;
0, 1./Ca;
0, 0];
M3 = [r(x2(t)-x1(t))./Rm; r(x1(t)-x4(t))./Ra];
dx = diff(x);
eqn = dx(t) == M1(t)*x(t) + M2*M3;
%now follow odeFunction first example
[eqs,vars] = reduceDifferentialOrder(eqn,x(t));
[M,F] = massMatrixForm(eqs,vars);
f = M\F;
odefun = odeFunction(f,vars);
%hen
initConditions = [vector of 5 values]; %[x1(0), x2(0), x3(0), x4(0), x5(0)];
tspan = [0 10]; %or as appropriate
mint = tspan(1);
maxt = tspan(end);
x0 = initConditions;
all_t = {};
all_x = {};
opts = odeset('events', @eventsFun);
while true
[t, x] = ode45(odefun, [mint, maxt], x0, opts);
all_t{end} = t;
all_x{end} = x;
mint = t(end);
x0 = x(end,:);
if mint == maxt; break; end
end
function [val, isterminal, dir] = eventsFun(t,x)
val(1) = x(2) - x(1);
val(2) = x(1) - x(4);
isterminal = [0; 0];
dir = [0; 0];
end
1 Comment
Jelena Kresoja
on 14 Aug 2020
This is the error that I get now:
Index exceeds matrix dimensions.
Error in odezero (line 60)
if (tL == t) && any(vL(indzc) == 0 & vR(indzc) ~= 0)
Error in ode45 (line 353)
odezero(@ntrp45,eventFcn,eventArgs,valt,t,y,tnew,ynew,t0,h,f,idxNonNegative);
Error in opet_jebeno (line 49)
[t, x] = ode45(odefun, [mint, maxt], x0, opts);
See Also
Categories
Find more on Calculus 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 (한국어)