Matlab ode15s change parameter value at specific time during solution
Show older comments
I'm trying to change the value of my variable Pin at specific points in time during the ode15s solution, in order to evaluate the dynamic response. But I get the error:
Error using odearguments (line 83)
The last entry in tspan must be different from the first entry.
I believe the error is somewhere here:
t_start=0;
t=t_start;
y=cond;
while idx_seg< length(t_seg)
idx_seg=idx_seg+1;
t_end=t_seg(idx_seg);
[t_sol,y_sol]=ode15s(@(t,y)f1_v1,[t_start, t_end],y(end,:));
t = [t; t_sol(2 : end)];
y = [y; y_sol(2 : end, :)];
t_start = t_end;
end
Here's the full code:
function [t,y]=f2_v4_21_08_18(cond, t_seg, Pin)
% -----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;
bp=0.8; ro=1040*10^6;
Ey=10^4*0.00750062;
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];
%%%%%
idx_seg=0;
function dy= f1_v1(t,y)
dy=zeros(13,1);
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);
% -----Constants -----
...
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];
for i=1:1:14
if i==1
pb(i)=Pin(idx_seg);
else
pb(i)=pb(i-1)-pdrop(i);
end
pb(i)=pb;
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
kin=v;
compl(j)=((3*pi*L(j)*r(j)^3)/(2*Ey*h(j)) )*10^-3;
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
%%differential eq
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];
dy=[dpt1;dpt2;dpt3;dpt4;dpt5;dpt6;dpt7;dpt8;dpt9;dpt10;dpt11;dpt12;dpt13];
end
t_start=0;
t=t_start;
y=cond;
while idx_seg< length(t_seg)
idx_seg=idx_seg+1;
t_end=t_seg(idx_seg);
[t_sol,y_sol]=ode15s(@(t,y)f1_v1,[t_start, t_end],y(end,:));
t = [t; t_sol(2 : end)];
y = [y; y_sol(2 : end, :)];
t_start = t_end;
end
for i=1:1:14
if i==1
pb=Pin(idx_seg);
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
kin=v;
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
for m=1:1:13
pt_weight(m)=y_sol(m)*Vt(m);
end
Vttot=Vt.*n;
Vt_sum=sum(Vttot);
ptot=sum(1/Vt_sum * (pt_weight));
end
These are the initial conditions and times I run it with:
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] = f2_v4_21_08_18(cond, [0, 10, 100], [60, 70, 60]);
plot(t, y);
Thanks in advance!
Accepted Answer
Torsten
on 22 Aug 2018
0 votes
t_start = 0 und t_seg(1) = 0.
In the first call to ODE15S, you consequently try to integrate from t_start=0 to t_end=0 which is not accepted by the integrator.
Best wishes
Torsten.
More Answers (0)
Categories
Find more on Programming in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
