Matlab ode15s在求解期间的特定时间更改参数值 [英] Matlab ode15s change parameter value at specific time during solution
问题描述
我正在尝试在ode15s解决方案的特定时间点更改变量Pin的值，以便评估动态响应.但是我得到了错误:
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.
我相信错误在这里:
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
这是完整的代码:
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);
提前谢谢！
推荐答案
检查while循环的第一次迭代:
Check the first iteration of your while loop:
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,:));
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,:));
会发生什么? idx_seg初始化为零，设置为1，即t_end设置为t_seg(1)(其第一个元素).您在调用f2_v4_21_08_18时将向量[0,10,100]传递给t_seg.因此，在您的第一次迭代中，t_end等于零.
What happens? idx_seg was initialized to zero, gets set to one, i.e., t_end is set to t_seg(1), its first element. You passed the vector [0,10,100] for t_seg in your call to f2_v4_21_08_18. Therefore, in your first iteration, t_end is equal to zero.
现在，t_start也是零，即您要ode15s运行从t_start到t_end的时间间隔，它们都是零.因此，它抱怨The last entry in tspan must be different from the first entry.
Now, t_start is also zero, i.e., you are asking ode15s to run the time interval from t_start to t_end which are both zero. Therefore, it complains that The last entry in tspan must be different from the first entry.
解决方案:传递不同的时间跨度向量，即[1, 10, 100]应该已经解决了该问题.您从使用[4,6,12]复制的示例.使它符合您的问题.
Solution: Pass a different time span vector, i.e., [1, 10, 100] should already solve the problem. The example you copied from used [4,6,12]. Make it match your problem.
编辑:再次检查了您在做什么，这是另一条建议:解释t_seg的方法是它包含参数更改的时间点.对于t_seg = [4,6,12]，这意味着从0到4，参数u(1)处于活动状态，从4到6，它是u(2)，从6到12，它是u(3).传递[0,10,100]和[60,70,60]意味着它将立即从60变为70，然后又返回60.这与传递[10,100]和[70,60]相同.
edit: Having checked again what you're doing, here is another piece of advice: the way to interpret t_seg is that it contains the time points where the parameter changes. For t_seg = [4,6,12] it means that from 0 to 4 the parameter u(1) is active, from 4 to 6 it is u(2) and from 6 to 12 it is u(3). Passing [0,10,100] and [60,70,60] would mean it changes from 60 to 70 right away, and then back to 60. This would be the same as passing [10,100] and [70,60].
这篇关于Matlab ode15s在求解期间的特定时间更改参数值的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
