When solve ODEs return NaN value
显示 更早的评论
Hallo,
I want to solve two sets of ODEs and use event function to swtich between them. And I also wish the parmeter of the coefficient matrix( crss and coss) of the ODEs will also change according to the results of ODEs in each step with Interp1();
In the first ODE seem work well. But when it switchs to the second, the ODEs cannot return values.
Can anyone helps me?
Thanks a lot.
clc;clear;
y0=[0; 0 ;0; 400; 0.9];
x=y0';
tspan=[0 1e-7];
tstart = tspan(1);
t=tstart;
tend=tspan(end);
vds1=linspace(0,400,1000);
crss=[];
coss=[];
for vds0=vds1;
if 0<=vds0 & vds0<=35
coss=[coss,1450e-12./((1+vds0/5).^0.8)];
elseif 35<vds0 & vds0<=400
coss=[coss,1450e-12./((1+vds0/1.3).^0.5)];
end
end
vds2=linspace(0,400,1000);
for vds0=vds2;
if 0<=vds0 & vds0<=20
crss=[crss,400e-12./((1+vds0/300).^40)];
elseif 20<vds0 & vds0<=35
crss=[crss,400e-12./((1+vds0/0.8).^0.8)];
elseif 35<vds0 & vds0<=400
crss=[crss,400e-12./((1+vds0/0.09).^0.49)];
end
end
ateout = [];
ayeout = [];
aieout = [];
fcn=@(t,x) odefun1(t,x,coss,vds1,crss,vds2);
options=odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@odefun1event);
while (t(end) < tend)
[at, ay, ate, aye, aie] = ode45(fcn, [t(end), tend], y0,options);
t = [t; at(2:end)];
x = [x; ay(2:end,:)];
y0 = x(end,:);
ateout = [ateout; ate];
ayeout = [ayeout; aye];
aieout = [aieout; aie];
if x(end,2) < 4
fcn=@(t,x) odefun1(t,x,coss,vds1,crss,vds2);
options=odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@odefun1event);
elseif x(end,3) < 10
fcn=@(t,x) odefun2(t,x,coss,vds1,crss,vds2);
options=odeset('AbsTol',1e-10,'RelTol',1e-12,'Events',@odefun2event);
else
break
end
end
plot(t,x(:,2),'-o')
function dx=odefun1(t,x,coss,vds1,crss,vds2)
dx=zeros(5,1);
Rg=5;
Lg=5e-9;
Ls=20e-9;
Ciss=1020e-12;
Coss=interp1(vds1,coss,x(4));
Crss=interp1(vds2,crss,x(4));
Cgs=Ciss-Crss;
V=x(4);
Vgs=15;
A=[-Rg/(Lg+Ls) -1/(Lg+Ls) 0 0 0; 1/Cgs 0 0 0 0; 0 0 0 0 0;0 0 0 0 0;0 0 0 0 0];
B=[Vgs/(Lg+Ls); 0;0;0;0];
dx=A*x+B;
end
function dx=odefun2(t,x,coss,vds1,crss,vds2)
dx=zeros(5,1);
Rg=5;
Lg=5e-9;
Ls=20e-9;
Coss=1000e-12;
Cgd=interp1(vds2,crss,x(4));
Ciss=1020e-12;
Cgs=Ciss-Cgd;
Vgs=15;
Vth=4;
Lp=100e-9;
Vf=0.9;
Rf=50e-3;
Vdc=400;
gf=10;
I0=10;
a=Lg+Ls-Ls^2/Lp;
b=Ciss*(Coss*Ciss-Cgd^2);
A=[-Rg/a -1/a Ls*Rf/(Lp*a) Ls/(Lp*a) 0;Coss*Ciss/b -gf*Cgd*Ciss/b Cgd*Ciss/b 0 0; Ls*Rg/(Lp*a) Ls/(Lp*a) -(Lg+Ls)*Rf/(Lp*a) -(Lg+Ls)/(Lp*a) 0; Cgd*Ciss/b -gf*Ciss^2/b Ciss^2/b 0 0;0 0 0 0 0];
B=[(Lp*Vgs-Ls*(Vdc+Rf*I0+Vf))/(Lp*a);gf*Cgd*Ciss*Vth/b; (-Ls*Vgs+(Lg+Ls)*(Vdc+Rf*I0+Vf))/(Lp*a);gf*Ciss^2*Vth/b;0];
dx=A*x+B;
end
function [value, isterminal, direction]= odefun1event(~,x)
value=x(2)-4;
isterminal=1;
direction=1;
end
function [value, isterminal, direction]= odefun2event(~,x)
value=x(3)-10;
isterminal=1;
direction=1;
end
4 个评论
Walter Roberson
2022-6-1
Using interp1 with the default interpolation method violates the continuity requirements of the mathematics behind Runge Kutta methods. You would need to use 'spline' or 'cubic' to be mathematically consistent.
Torsten
2022-6-1
Having in mind how ODE-solvers are misused in this forum: this is a harmless one.
Walter Roberson
2022-6-1
When the context is continuous control and the interval between force adjustments are not so far apart and the adjustments are relatively smooth, the result can be sometimes be "good enough". A real car might take a corner a bit more sharply or slow a bit sooner, and that might not matter. But if you are studying pumping brakes, or you have geers with fine tolerance, the error might be too much.
Torsten
2022-6-1
Linear interpolation is even taken as a MATLAB example on how to deal with time-dependent terms in ODEs:
Example "ODE with Time-Dependent Terms" under
Should be deleted -:)
采纳的回答
The command
Cgd=interp1(vds2,crss,x(4));
in odefun2 returns NaN.
The reason is that x(4) > 400 and you only have specified crss for vds2 <= 400. This means interp1 has to extrapolate outside the interval where values for crss have been specified. Since you didn't specify "extrap" in the call to interp1, the default method is chosen and NaN is returned.
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
产品
Translated by 
