2nd order differential eqn for Windkessel model

12 views (last 30 days)
I am struggling to solve this 2nd order ODE in matlab for a 4WK Parallel Windkessel model, this is what I have so far for my input and outputs for a 2WK model, a 3WK model (attached).
I am trying to solve a 2nd order ODE:
I''(R*L*C*r) + I'(L*(R+r)) +I*(R*r) = Y"(L*C*R)+Y'(C*R*r +L)+ Y*r
I is :
I=@(t)I0*sin((pi*t)/Ts).^2.*(t<=Ts); %input current flow
t=0:0.001:Tc;
Initial conditions for Y(0)= 80, the rest of the variables are known constants.
Thanks in advance!
  7 Comments
Nang Su Lwin
Nang Su Lwin on 22 Nov 2022
Edited: Nang Su Lwin on 22 Nov 2022
I tried to simplify the equations as much as possible to help me not go crazy as I tried to solve this equation (I substituted in all the constants so I didn't have to deal with so many variables). The problem is that I don't know how to implement the fact that my Input function I,I', and I'' (now to the right of the equation) doesn't consider that they should only output when t<0.33333:
syms y(t)
Dy = diff(y,t);
D2u = diff(y,t,2);
ode= 0.05*y+0.055*diff(y,t,1)+ +0.005*diff(y,t,2) == 15.75*pi*cos(3*pi*t)*sin(3*pi*t)+22.5*pi^2*cos(3*pi*t)^2 - 22.5*pi^2*sin(3*pi*t)^2 +25*sin(3*pi*t).^2;
cond1 = y(0) == 80;
cond2 = Dy(0) == 0;
conds = [cond1 cond2];
ySol(t) = dsolve(ode, conds);
ySol= simplify(ySol)
t=0:0.001:60/72;
plot(t,ySol(t),'g')
So I get something like this:
Instead of what I'm expecting:

Sign in to comment.

Accepted Answer

Torsten
Torsten on 22 Nov 2022
Edited: Torsten on 22 Nov 2022
I = @(t) I0*sin((pi*t)/Ts).^2.*(t<=Ts); %input current flow
Idot = @(t) I0*2*sin(pi*t/Ts).*cos(pi*t/Ts)*pi/Ts.*(t<=Ts);
Idotdot = @(t) I0*2*(cos(pi*t/Ts).^2 - sin(pi*t/Ts).^2)*(pi/Ts)^2.*(t<=Ts);
fun = @(t,y)[y(2);(Idotdot(t)*(R*L*C*r)+Idot(t)*(L*(R+r))+I(t)*(R*r) - (y(2)*(C*R*r+L)+y(1)*r))/(L*C*R)];
y0 = [80 ;0];
tspan = [0 5/6];
[T,Y] = ode45(fun,tspan,y0)
  2 Comments
Nang Su Lwin
Nang Su Lwin on 22 Nov 2022
Edited: Torsten on 22 Nov 2022
Just for others who might read this post later and are trying to also plot some Windkessel responses, this is the 4Windkessel code:
clear all
clc
I0= 500; % maximum flow
Tc=60/72; % heart period
Ts=(2/5*Tc); % time in systole
P_ss=80; % diastolic pressure
R= 1;
C=1;
R1=0.05;
L=0.005;
I=@(t)I0*sin((pi*t)/Ts).^2.*(t<=Ts); %input current flow
Idot = @(t)I0*2*sin(pi*t/Ts).*cos(pi*t/Ts)*pi/Ts.*(t<=Ts);
Idotdot = @(t) I0*2*(cos(pi*t/Ts).^2 - sin(pi*t/Ts).^2)*(pi/Ts)^2.*(t<=Ts);
fun = @(t,y)[y(2);(Idotdot(t)*(R*L*C*R1)+Idot(t)*(L*(R+R1))+I(t)*(R*R1) - (y(2)*(C*R*R1+L)+y(1)*R1))/(L*C*R)];
y0 = [80 ;0];
tspan = [0 60/72];
[T,Y] = ode45(fun,tspan,y0);
plot(T,Y(:,1),'g')

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!