Is it possible to get a time-dependent function as an output of an ODE solver?

조회 수: 1 (최근 30일)
Nang Su Lwin
Nang Su Lwin 2022년 11월 23일
답변: Steven Lord 2022년 11월 24일
My code is below, it plots as I expect it to. I would like to know if it is possible to get an output where the ODE solver gives the solution as a function in time rather than 2 column vectors [T,Y]. I wanted to use the output to further analyze something, but need a function in time rather than the column vectors.
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')

답변 (2개)

Torsten
Torsten 2022년 11월 23일
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')
function_of_time = @(t)interp1(T,Y(:,1),t);
function_of_time(30/72)
ans = 119.2901

Steven Lord
Steven Lord 2022년 11월 24일
Call ode45 with one output argument. Pass that output from ode45 and an array of times into the deval function to evaluate the solution at those times.

카테고리

Help CenterFile Exchange에서 Ordinary Differential Equations에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by