Matlab ode15s change parameter value at specific time during solution
1 回表示 (過去 30 日間)
古いコメントを表示
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!
0 件のコメント
採用された回答
Torsten
2018 年 8 月 22 日
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.
2 件のコメント
Torsten
2018 年 8 月 22 日
cond is 13x1, thus y is 13x1 at the start. Consequently y(end,:) is a scalar, not a vector of length 13.
その他の回答 (0 件)
参考
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!