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.

[Walter Roberson]

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]

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]

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]

So, C=1/E, and Cdot I also have to compute. In the eqution below tn=t/Tmax where Tmax is know.

[Walter Roberson]

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 ??

[Jelena Kresoja]

E is supposed to be a vector.

[Sam Chak]

Please confirm if your En(tn) looks as such:

What are your Emin and Emax?

[Jelena Kresoja]

Yes, this is it. Emax=2, Emin=0.6

[Jelena Kresoja]

Typo, Emin=0.06, sorry

[Walter Roberson]

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]

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]

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]

You can follow the advices from Walter Roberson and James Tursa, or refer to the examples in this link:

https://www.mathworks.com/help/matlab/ref/ode45.html

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]

Thank you! I will try to solve it!

[Jelena Kresoja]

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]

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]

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]

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]

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]

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]

Ra and Rm are just constants. The error is for M2*M3. The size of 1x1 and for M3 its 10x1.

[Walter Roberson]

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]

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.

[Walter Roberson]

post your code attempt.

[Jelena Kresoja]

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]

please format as code. I cannot copy and paste that from my mobile phone.

[Jelena Kresoja]

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]

 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]

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.

[Jelena Kresoja]

Its still not working unfortunatelly. Still the same error message.

[Walter Roberson]

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]

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]

You must use event functions for that situation, as it is not differentiateable at the boundary.

[Sam Chak]

Hi Jelena Kresoja

Say . The desired piecewise function as described is given by

Good luck!

[Walter Roberson]

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]

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?

[Walter Roberson]

Not even close. Please review https://www.mathworks.com/help/matlab/math/ode-event-location.html

[Jelena Kresoja]

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]

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]

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]

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]

Im so sorry i literally dont get whats written here, or where im supposed to put this.....

Answer 1

[hosein Javan]

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.

Answer 2

[Walter Roberson]

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

[Jelena Kresoja]

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);