Implicit system of ODEs, problem with ode15i

5 次查看（过去 30 天）
Stephane
Stephane 2023-5-13
评论： Stephane 2023-5-15
Hello,
I am trying to solve a system of two 2nd-order ODEs that I write as four 1st-order ODEs. The equations are quite non-linear and non-autonomous (there is a periodic forcing).
I don't know if my code is useful but I am including it below. Theres quite a few parameters, I am only changing the one I call Gamma in my code. This is my main script:
clc; clear all; close all
%% constant parameters
global B Gs Gamma I k2 C100 C001 C011 C021 C031 f2_1
k2 = 2.36502037243135201301202405;
C100 = 0.20253735124738893522296764; C001 = -2.2937835087802942891067795;
C011 = 0.6416819172490610458168773; C021 = 1.185226949697087841616967;
C031 = 0.894260611316773729819612; f2_1 = -1.4266400021183172634428127;
B=1e-4; I=1e-2; Gs=1e-4;
Gamma=60;
%% solving
t0=0; tend=200;
tspan=0:0.01:tend;
options = odeset('RelTol',1e-6,'AbsTol',1e-6);
% wrong initial guess
y0 = [1 0 0 0]; yp0 = [0 0 0 0];
% computes an initial guess
[y0,yp0] = decic(@odefun,t0,y0,[1 0 0 0],yp0,[],options);
[t,y] = ode15i(@odefun,tspan,y0,yp0);
a0=y(:,1); a0p=y(:,2);
a2=y(:,3); a2p=y(:,4);
%% plotting
figure();
plot(t,a0,'-o')
title('a0, y(:,1)')
xlabel('t'); ylabel('y(:,1)')
figure();
plot(t,a2)
title('a2, ty(:,2)')
xlabel('t'); ylabel('y(:,3)')
The ODEs are define by the following, where the global variables are simply constant parameters:
function res = odefun(t,y,yp)
global B Gs Gamma I k2 C100 C001 C011 C021 C031 f2_1
res = zeros(4,1);
res(1) = yp(1)-y(2);
res(2) = yp(3)-y(4);
res(3) = I*yp(2)+(B*k2^4*y(3)+I*yp(4))*f2_1 + Gs + Gs*Gamma*cos(t);
res(4) = -0.5*y(2)-C100*y(4) - ((I*yp(2)-Gs-Gs*Gamma*cos(t))/f2_1) .* ...
( C001*y(1).^3 + 3*C011*y(1).^2.*y(3) + 3*C021*y(1).*y(3).^2 + C031*y(3).^3 );
end
The ODEs have a forcing term in cos(t) that is proportional to the parameter Gamma.
For Gamma<=60, the solution looks like what I would expect physically (image below is Gamma=60):
However upon increasing Gamma, it looks like some sort of singularity appears. Putting Gamma=70, we can see a sharp change in the solution near t=15:
For Gamma=80 something similar occurs near t=2.8, and the solver actually completely fails early:
I suspect that this singularity for high values of Gamma is not a feature of the equations but stems from numerical errors. How could I investigate this more? This is my first time solving implicit ODEs and there does not seem to be many alternatives to ODE15i. Changing the solver tolerance doesn't really do much, and I don't know what else to try.
EDIT: It appears that the initial condition is criticial, and minute changes matter. For instance with Gamma=70, adding the following line after calling the decic function completely changes the solution:
yp0(4)=yp0(4)+1e-15;
I overlooked this before and I actually dont really know what decic does ; I'll read a bit on this.

请先登录，再回答此问题。

回答（1 个）

Torsten
Torsten 2023-5-13
编辑：Torsten 2023-5-13
You can explicitly solve for yp - so there is no need to use ode15i. But the Gamma-problem remains.
clc; clear all; close all
%% constant parameters
global B Gs Gamma I k2 C100 C001 C011 C021 C031 f2_1
k2 = 2.36502037243135201301202405;
C100 = 0.20253735124738893522296764; C001 = -2.2937835087802942891067795;
C011 = 0.6416819172490610458168773; C021 = 1.185226949697087841616967;
C031 = 0.894260611316773729819612; f2_1 = -1.4266400021183172634428127;
B=1e-4; I=1e-2; Gs=1e-4;
Gamma=60;
%% solving
t0=0; tend=200;
tspan=0:0.01:tend;
options = odeset('RelTol',1e-8,'AbsTol',1e-8);
% wrong initial guess
y0 = [1 0 0 0];
% computes an initial guess
[t,y] = ode15s(@odefun,tspan,y0);%,options);
a0=y(:,1); a0p=y(:,2);
a2=y(:,3); a2p=y(:,4);
%% plotting
figure();
plot(t,a0,'-o')
title('a0, y(:,1)')
xlabel('t'); ylabel('y(:,1)')
figure();
plot(t,a2)
title('a2, ty(:,2)')
xlabel('t'); ylabel('y(:,3)')
function yp = odefun(t,y)
global B Gs Gamma I k2 C100 C001 C011 C021 C031 f2_1
yp = zeros(4,1);
yp(1) = y(2);
yp(3) = y(4);
yp(2) = ((-0.5*y(2)-C100*y(4))*f2_1/...
( C001*y(1).^3 + 3*C011*y(1).^2.*y(3) + 3*C021*y(1).*y(3).^2 + C031*y(3).^3 )+...
(Gs+Gs*Gamma*cos(t)))/I;
yp(4) = ((-I*yp(2)-Gs-Gs*Gamma*cos(t))/f2_1 - B*k2^4*y(3))/I;
% res(1) = yp(1)-y(2);
% res(2) = yp(3)-y(4);
% res(3) = I*yp(2)+(B*k2^4*y(3)+I*yp(4))*f2_1 + Gs + Gs*Gamma*cos(t);
% res(4) = -0.5*y(2)-C100*y(4) - ((I*yp(2)-Gs-Gs*Gamma*cos(t))/f2_1) .* ...
% ( C001*y(1).^3 + 3*C011*y(1).^2.*y(3) + 3*C021*y(1).*y(3).^2 + C031*y(3).^3 );
end
  1 个评论
Stephane
Stephane 2023-5-15
This is true, I'm not sure why I got confused and thought I needed ode15i... Thanks! It will be helpful to be able to use more standard ode solvers to try to understand what's going on

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

标签

产品


版本

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by