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

[Walter Roberson]

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]

Having in mind how ODE-solvers are misused in this forum: this is a harmless one.

[Walter Roberson]

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]

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

https://de.mathworks.com/help/matlab/ref/ode45.html

Should be deleted -:)

Answer 1

[Torsten]

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.

[Tianyu Li]

Hi Torsten,

I have fix it with your answers. Thanks a lot.