ode15i instability and initial values

1 view (last 30 days)
Jan Filip
Jan Filip on 18 Feb 2021
Edited: Jan Filip on 18 Feb 2021
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?

Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!