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 次查看(过去 30 天)
显示 更早的评论
42 个评论
Walter Roberson
2020-8-4
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
2020-8-4
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
2020-8-4
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
2020-8-5
So, C=1/E, and Cdot I also have to compute. In the eqution below tn=t/Tmax where Tmax is know.
Walter Roberson
2020-8-5
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
2020-8-5
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
2020-8-6
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
2020-8-6
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
2020-8-6
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
2020-8-8
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
2020-8-8
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
2020-8-9
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
2020-8-9
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
2020-8-10
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
2020-8-10
编辑:Walter Roberson
2020-8-10
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
2020-8-10
Ra and Rm are just constants. The error is for M2*M3. The size of 1x1 and for M3 its 10x1.
Walter Roberson
2020-8-10
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
2020-8-11
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
2020-8-11
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;
Jelena Kresoja
2020-8-11
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
2020-8-11
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
2020-8-11
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
2020-8-11
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
2020-8-11
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
2020-8-11
You must use event functions for that situation, as it is not differentiateable at the boundary.
Walter Roberson
2020-8-12
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
2020-8-12
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
2020-8-13
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
2020-8-13
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
2020-8-14
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
2020-8-14
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
2020-8-14
Im so sorry i literally dont get whats written here, or where im supposed to put this.....
回答(2 个)
hosein Javan
2020-8-10
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
2020-8-14
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 个评论
Jelena Kresoja
2020-8-14
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);
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!发生错误
由于页面发生更改,无法完成操作。请重新加载页面以查看其更新后的状态。
您也可以从以下列表中选择网站:
如何获得最佳网站性能
选择中国网站(中文或英文)以获得最佳网站性能。其他 MathWorks 国家/地区网站并未针对您所在位置的访问进行优化。
美洲
- América Latina (Español)
- Canada (English)
- United States (English)
欧洲
- 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)
亚太
- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)