Matlab ode15s change parameter value at specific time during solution

조회 수: 6(최근 30일)
gorilla3
gorilla3 2018년 8월 22일
댓글: Torsten 2018년 8월 22일
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!

채택된 답변

Torsten
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
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!

Translated by