ode15i instability and initial values
1 view (last 30 days)
Show older comments
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?
0 Comments
Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!