MATLAB: Attempting to solve differential equation model of windkessel two element using ode45, getting errors

ode 45

Hello, I am trying to use ode45 to solve numerically for pressure for the equation below: dP/dt = Q(t)/C -P(t)/R*C where R and C are known constants and I am given Q(t) as a matrix of size (206X1) stored in variable named flow. I'm suppose to use Q(t) to solve numerically for pressure and I tried using ode45 but I get the error "Error using odearguments (line 92) @(T,Y)QT/C-Y/R*C returns a vector of length 206, but the length of initial conditions vector is 1. The vector returned by @(T,Y)QT/C-Y/R*C and the initial conditions vector must have the same number of elements." this is my code: tspan= 0:1:206 Qt = flow(:,:) y0 = 0; [t,y] = ode45(@(t,y) Qt/C-y/R*C,tspan,y0) It seems my problem is with figuring out how to put Q(t) in the write syntax for ode45.

Thanks.

Best Answer

You are doing curve fitting (also known as parameter estimation). The code in Monod kinetics and curve fitting (link), although an enzyme kinetics problem, will demonstrate how to integrate a differential equation using ode45 and estimate the parameters of the system. You should easily be able to adapt it to your problem.

Related Solutions

MATLAB: System identification Windkessel – pulse wave

If you provide the input and output as a ‘.mat’ file, and your model, it should be possible to derive the component values. (I can derive a circuit model for it if you want.)

You aren’t doing system identification since you already have a model for your system. You’re doing parameter estimation of your model.

EDIT —

The circuit model I derived for your model (I let the Symbolic Math Toolbox do the algebra) is:

syms C L R r Pi Po P2 s t 
Z1 = 1/r;
Z2 = s*C + 1/R;
Z3 = 1/(s*L);
Zt = Z1 + Z2 + Z3;
I = Pi/Zt;
Po = Z2*I;
Out = ilaplace(partfrac(Po, s));
Out = simplify(Out, 'Steps',20);
Out_fcn = matlabFunction(Out, 'Vars',{t,Pi,r,C,R,L})
Out_fcn = @(t,Pi,r,C,R,L) Pi.*dirac(t)-(Pi.*exp((t.*(R+r).*(-1.0./2.0))./(C.*R.*r)).*(cosh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r))-1.0./sqrt(L).*sinh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r)).*(L.*r+L.*R-C.*R.*r.^2.*2.0).*1.0./sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0)))./(C.*r);

The ‘Z’ values are the three lumped impedances (and ‘Zt’, the total impedance), where the output is taken across ‘Z2’. The idea is straightforward: the current ‘I’ (or blood flow here) through the series of impedances is simply ‘Pi/Zt’. It is then straightforward to calculate the voltage (or pressure) drop across ‘Z2’ as ‘Z2*I’.

In order to fit your function and estimate the parameters, we need three data vectors: ‘t’, ‘Pi’, and ‘Po’, the time, input pressure with respect to time, and output pressure with respect to time. The function fits the output pressure it estimates with specific parameters to the ‘Po’ vector.

The procedure for estimating the parameters would then be:

% % % MAPPING: b(1) = r,  b(2) = C,  b(3) = R,  b(4) = L
Out_fcn = @(t,Pi,r,C,R,L) Pi.*dirac(t)-(Pi.*exp((t.*(R+r).*(-1.0./2.0))./(C.*R.*r)).*(cosh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r))-1.0./sqrt(L).*sinh((1.0./sqrt(L).*t.*sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0).*(1.0./2.0))./(C.*R.*r)).*(L.*r+L.*R-C.*R.*r.^2.*2.0).*1.0./sqrt(L.*R.^2+L.*r.^2-C.*R.^2.*r.^2.*4.0+L.*R.*r.*2.0)))./(C.*r);
t = t(:);                                                           % Column Vector


Pi = Pi(:);                                                         % Column Vector
Po = Po(:);                                                         % Column Vector
b0 = ['Initial Parameter Estimates For r,C,R,L — IN THAT ORDER — As A Vector'];
nrmrsd = @(b) norm(Po - Out_fcn(t,Pi,b(1),b(2),b(3),b(4)));         % Norm Of Residuals (Errors)
bvals = fminsearch(nrmrsd, b0);

The ‘bvals’ vector will be the estimated parameters.

NOTE — This is UNTESTED CODE since I do not have your data to test it with.

MATLAB: Plot just y or y’ in a second order eq.

syms y(t) x(t)
[V] = odeToVectorField(diff(y, 2) == -1*(1-1*cos(2*t))*y);
M = matlabFunction(V,'vars', {'t','Y'});
[ty,y] = ode45(@(t,Y)M(t,Y),tspan,[0 1]);
[tx,x] = ode45(@(t,Y)M(t,Y),tspan,[1 0]);
hold on
plot(tx, x(:,1));
plot(ty, y(:,2));
hold off
legend({"x", "y'"})

Best Answer

Related Solutions

MATLAB: System identification Windkessel – pulse wave

MATLAB: Plot just y or y’ in a second order eq.

Related Question