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.

3 次查看（过去 30 天）
Jelena Kresoja
Jelena Kresoja 2020-8-4
评论： Jelena Kresoja 2020-8-14
  42 个评论
显示 40更早的评论
Walter Roberson
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()

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

hosein Javan
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
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
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!

Translated by