Implicit system of ODEs, problem with ode15i

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.

Answer 1

[Torsten]

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

[Stephane]

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