ode15i instability and initial values

1 回表示 (過去 30 日間)
jf
jf 2021 年 2 月 18 日
編集済み: jf 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?

回答 (0 件)

カテゴリ

Help Center および File ExchangeMathematics についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by