ode15s extract solutions for every time step

조회 수: 15(최근 30일)
gorilla3
gorilla3 2018년 8월 23일
답변: Steven Lord 2018년 8월 23일
Hi,
I'm running this system of ode with ode15s, and I would like to visualize one of the parameters ( ptot) at each step of the solution ( t). The timespan of the ode is 0:0.1:100, so I should have 1001 points. However, the .txt files where I'm printing ptot and t only have 44 elements each. Could you please help me find the mistake? I'm trying to plot ptot against t because I'm changing a parameter ( Pin) during the solution and would like to see the effect of this change in time.
Here's the code:
tspan=0:0.1:100;
cond=[51.2112 ; 63.8766 ; 60.7979 ; 49.0010 ; 35.3767 ; 28.5718 ; 33.7930 ; 31.1300 ; 30.6594 ; 29.9741 ; 30.2541 ; 29.6828 ; 28.9798 ];
[t,y, ptot] = ode15s(@fun2_23_08_18_v2,tspan,cond,[])
Ptot = 'ptot_ch.txt';
ttot = 't_ch.txt';
pt_read = importdata(Ptot);
tread=importdata(ttot);
plot(tread,pt_read);
%%% FUNCTION
function [dz,ptot]= fun2_23_08_18_v2(t,y)
pt1=y(1); pt2=y(2); pt3=y(3); pt4=y(4); pt5=y(5); pt6=y(6); pt7=y(7); pt8=y(8); pt9=y(9); pt10=y(10); pt11=y(11); pt12=y(12); pt13=y(13);
if t<10
Pin=60;
else
Pin=70;
end
% -----Constants -----
N=3.38*10^6; k=2.96*10^-7; alphat=2.6*10^-5; chb=0.2; M=30*10^-9; K=5*10^(-8)*10^-3; H=0.42; S0=0.98;
v=3* 7.5*10^-6;
r=[0.0119850000000000;0.00958500000000000;0.00764000000000000;0.00604000000000000;0.00473000000000000;0.00366000000000000;0.00400000000000000;0.00575500000000000;0.00726500000000000;0.00889500000000000;0.0107250000000000;0.0128500000000000;0.0153850000000000]; %mm
L=[1.27076497943190;0.932650928622621;0.544932761536915;0.303082765473283;0.161799106136796;0.155424891414508;0.245072221621871;0.475103125625241;0.273016623935407;0.427646038844292;0.634082325832342;0.846354695529459;0.938696601022114]; %mm
h=[0.00484000000000000;0.00425000000000000;0.00381000000000000;0.00349000000000000;0.00327000000000000;0.00314000000000000;0.000309000000000000;0.00115000000000000;0.00145000000000000;0.00178000000000000;0.00215000000000000;0.00257000000000000;0.00308000000000000];%mm
mu=[1.19409872390289e-05;1.12760032450214e-05;1.06583134073916e-05;1.00804835896938e-05;9.56162410894012e-06;9.20633512007761e-06;9.29628357913371e-06;9.96996247072375e-06;1.05291656347798e-05;1.10660983739492e-05;1.16035344790804e-05;1.21594980256614e-05;1.27473361251949e-05]; %mmHg*s
R= [1.8728 3.1728 4.3411 5.8457 7.8705 20.3059 22.6623 10.9961 2.6277 1.9250 1.4161 0.9612 0.5439]; %[mmHg*s/ml]
pdrop=[0,6.93,5.87,4.02,2.70,1.82,2.35,2.62,1.27,0.61,0.89,1.31,1.78,2.01];
n=[1;2;4;8;16;32;64;32;16;8;4;2;1];
for i=1:1:14
if i==1
pb(i)=Pin;
else
pb(i)=pb(i-1)-pdrop(i);
end
end
diffp=diff(pb)*(-1);
pb_s=[pb(1)+pb(2);pb(2)+pb(3);pb(3)+pb(4);pb(4)+pb(5);pb(5)+pb(6);pb(6)+pb(7);pb(7)+pb(8);pb(8)+pb(9); pb(9)+pb(10);pb(10)+pb(11);pb(11)+pb(12);pb(12)+pb(13);pb(13)+pb(14)];
for z=1:1:14
S(z)=(((pb(z)^3+150*pb(z))^(-1) *23400)+1)^(-1);
end
deltaS=diff(S)*(-1);
for j=1:1:13
compl(j)=((3*pi*L(j)*r(j)^3)/(2*Ey*h(j)) )*10^-3; %vessel compliance (ElBouri2018) [ml/mmHg]
Vb(j)=(compl(j)/2)*pb_s(j); %[ml]
q(j)=diffp(j)/R(j);
Vt(j)= chb*H*q(j)*deltaS(j)/M;
end
dpt1=1/(alphat*Vt(1))*( (2*pi*K*r(1)*L(1))/h(1) *(1/2*(pb(1)+pb(2)) -pt1) -M*Vt(1));
dpt2=1/(alphat*Vt(2))*( (2*pi*K*r(2)*L(2))/h(2) *(1/2*(pb(2)+pb(3)) -pt2) -M*Vt(2));
dpt3=1/(alphat*Vt(3))*( (2*pi*K*r(3)*L(3))/h(3) *(1/2*(pb(3)+pb(4)) -pt3) -M*Vt(3));
dpt4=1/(alphat*Vt(4))*( (2*pi*K*r(4)*L(4))/h(4) *(1/2*(pb(4)+pb(5)) -pt4) -M*Vt(4));
dpt5=1/(alphat*Vt(5))*( (2*pi*K*r(5)*L(5))/h(5) *(1/2*(pb(5)+pb(6)) -pt5) -M*Vt(5));
dpt6=1/(alphat*Vt(6))*( (2*pi*K*r(6)*L(6))/h(6) *(1/2*(pb(6)+pb(7)) -pt6) -M*Vt(6));
dpt7=1/(alphat*Vt(7))*( (2*pi*K*r(7)*L(7))/h(7) *(1/2*(pb(7)+pb(8)) -pt7) -M*Vt(7));
dpt8=1/(alphat*Vt(8))*( (2*pi*K*r(8)*L(8))/h(8) *(1/2*(pb(8)+pb(9)) -pt8) -M*Vt(8));
dpt9=1/(alphat*Vt(9))*( (2*pi*K*r(9)*L(9))/h(9) *(1/2*(pb(9)+pb(10)) -pt9) -M*Vt(9));
dpt10=1/(alphat*Vt(10))*( (2*pi*K*r(10)*L(10))/h(10) *(1/2*(pb(10)+pb(11)) -pt10) -M*Vt(10));
dpt11=1/(alphat*Vt(11))*( (2*pi*K*r(11)*L(11))/h(11) *(1/2*(pb(11)+pb(12)) -pt11) -M*Vt(11));
dpt12=1/(alphat*Vt(12))*( (2*pi*K*r(12)*L(12))/h(12) *(1/2*(pb(12)+pb(13)) -pt12) -M*Vt(12));
dpt13=1/(alphat*Vt(13))*( (2*pi*K*r(13)*L(13))/h(13) *(1/2*(pb(13)+pb(14)) -pt13) -M*Vt(13));
pt_tot=[pt1;pt2;pt3;pt4;pt5;pt6;pt7;pt8;pt9;pt10;pt11;pt12;pt13];
for l=1:1:13
pt_weight(l)=pt_tot(l)*Vt(l);
Vtt(l)=Vt(l)*n(l);
end
Vt_sum=sum(Vtt);
ptot=sum(1/Vt_sum * (pt_weight));
dz=[dpt1;dpt2;dpt3;dpt4;dpt5;dpt6;dpt7;dpt8;dpt9;dpt10;dpt11;dpt12;dpt13];
%STORE VALUES
fileID=fopen('ptot_ch.txt','a');
fprintf(fileID,' %2.9f\n',ptot);
fclose(fileID);
fileID=fopen('t_ch.txt','a');
fprintf(fileID,' %2.9f\n',t);
fclose(fileID);
fileID=fopen('Vt.txt','a');
fprintf(fileID,' %2.9f\n',Vt);
fclose(fileID);

채택된 답변

Steven Lord
Steven Lord 2018년 8월 23일
There's no mistake. From the description of the tspan input argument on the documentation page for the ode15s function:
"If tspan has more than two elements [t0,t1,t2,...,tf], then the solver returns the solution evaluated at the given points. However, the solver does not step precisely to each point specified in tspan. Instead, the solver uses its own internal steps to compute the solution, then evaluates the solution at the requested points in tspan. The solutions produced at the specified points are of the same order of accuracy as the solutions computed at each internal step."
To anticipate your next question no, there is no way to force ode15s to step exactly to the list of points in your tspan vector. If it calculates that it can satisfy your tolerances while taking larger (and fewer) steps than you specified it will.
Since you want the value of a variable that's not returned from your function to the ODE solver, the easiest way to do that is probably to modify your function to return that value as the third output argument, call ode15s normally, then call your function using the solution from ode15s as input while specifying the third output. This also avoids needing to worry about rejected steps, where the solver realizes it took too long a step and evaluates your function with a smaller step from the previous point. This appears like the solver stepped backwards, but instead it effectively threw away the longer step and replaced it with the smaller.

추가 답변(0개)

Community Treasure Hunt

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

Start Hunting!

Translated by