ode15i instability and initial values
3 次查看(过去 30 天)
显示 更早的评论
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/523462/image.png)
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 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Assembly 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!