ode15i instability and initial values

2 次查看（过去 30 天）

Jan Filip 2021-2-18

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/748867-ode15i-instability-and-initial-values

编辑： Jan Filip 2021-2-18

Consider following MWE:

eqn1 =  i_V11(t) + u_1(t)/10 - u_2(t)/10 == 0
eqn2 = (3*u_2(t))/25 - u_1(t)/10 - u_3(t)/50 + (7737125245533627*diff(u_2(t), t))/154742504910672534362390528 - (7737125245533627*diff(u_4(t), t))/154742504910672534362390528 == 0
eqn3 =                                                                                  u_3(t)/50 - u_2(t)/50 + (7737125245533627*diff(u_3(t), t))/154742504910672534362390528 == 0
eqn4 =                              i_L11(t) - (7737125245533627*diff(u_2(t), t))/154742504910672534362390528 + (7737125245533627*diff(u_4(t), t))/154742504910672534362390528 == 0
eqn5 =                                                                                                                                    diff(i_L11(t), t)/100000000 - u_4(t) == 0
eqn6 =                                                                                                                                                                  u_1(t) == (-0.012*sin(2*pi*1e6*t) * 0.5*sin(3*pi*1e6*t))                                                                                                                                                     
eqns = [eqn1; eqn2; eqn3; eqn4; eqn5; eqn6];
newEqs = reduceDAEToODE(eqns, x)
M = incidenceMatrix(eqns,x)
isLowIndexDAE(eqns,x)
[DAEs,DAEvars] = reduceDAEIndex(eqns,x)
f = daeFunction(DAEs,DAEvars);
F = @(t,Y,YP) f(t,Y,YP);
y0est = [0 0 0 0 0 0];
yp0est = zeros(size(DAEvars,1),1);
[y0,yp0] = decic(F,0,y0est,[],yp0est,[]);
%case1
tspan = [0 9e-6];
[tSol,ySol] = ode15i(F,tspan,y0,yp0);
figure
plot(tSol,ySol(:,1))
hold on
%case 2
tspan = [0 10e-6];
[tSol,ySol] = ode15i(F,tspan,y0,yp0,'r');
plot(tSol,ySol(:,1))
%case 3
tspan = [10e-9 100e-6];
[tSol,ySol] = ode15i(F,tspan,y0,yp0);
plot(tSol,ySol(:,1),'g:')
legend('tspan 0:9e-6','tspan 0:10e-6','tspan 10e-9:100e-6')
xlim([0 9e-6])

For a given DAE system wirth consistent intiial value

y0est = [0 0 0 0 0 0];

ode15i solver gives right answer on span [0 9e-6]. For longer span [0 10e-6] however it fails terribly. Far more confusing is the behaviour on span [10e-9 100e-6]. It is of course inconsistent, but solution is stable for large spans reaching [10e-9 1] and even more.

See figure below.

What is wrong with integration from 0 and how to avoid such behaviour?

EDIT: using explicitly given time points in tspan = [0:1e-9:400e-6] is useful. Integration then works as desired. Still, can someone explain this behaviour and possibly give some treatment tips?