Solving nonlinear system with varying parameter

Question

I have the Jansen-Rit system of 6 ode. I am trying to solve it for a range of values of one of its parameters, parameter 'p', using a for loop at the integration.

I have managed to solve this Jansen-Rit for all 'p' values using ode45, however, using a non-built-in integrator, the results are not correct. The rank of the calculated matrices with the 'integrator' is way lower (2) than the rank therein while using ODE45 (60), for identical initial conditions, parameter values and time-step.

Using this 'integrator' though, I have managed to solve correctly the Van der Pol system, for a range of values of its parameter.

I would definitely appreciate any feedback on what I am doing wrong

Using ODE45 (WORKS)

p = -100:+1:350;
time = 0:.05:3;
initial = [0 0 0 0 0 0];
x = NaN(length(time),length(initial),length(p));

for i=1:length(p)

    [t,x(:,:,i)] = ode45(@ode,time,initial,[],p(i));

end
figure(1)
plot3(squeeze(x(:,1,:)),squeeze(x(:,2,:)),squeeze(x(:,3,:)))

function dx = ode(~,x,p)

A = 3.25e-3;
B = 22e-3;
a = 100;
b = 50;
C = 135;
C1 = C;
C2 = 0.8*C;
C3 = 0.25*C;
C4 = 0.25*C;
dx = NaN(6,1);
dx(1) = x(4);
dx(2) = x(5);
dx(3) = x(6);
dx(4) = A*a*sigm(x(2)-x(3)) - 2*a*x(4) - x(1)*a^2;
dx(5) = A*a*(p+C2*sigm(C1*x(1))) - 2*a*x(5) - x(2)*a^2;
dx(6) = B*b*C4*sigm(C3*x(1)) - 2*b*x(6) - x(3)*b^2;
end

function X = sigm(u)
u0 = 6e-3;
e0 = 2.5;
r = 0.56e3;
X = 2*e0/(1+exp(r*(u0-u)));
end

Using the non-built-in integrator (DOES NOT WORK)

figure
hold all

for p = -100:+1:350
    inicond = [0 0 0 0 0 0];
    dt = 0.01;
    time = 0:dt:10;

    [y] = integrator(@ode,inicond,time,dt,p);
    x1 = y(:,1);
    x2 = y(:,2);
    x3 = y(:,3);
    figure(1)
    plot3(x1,x2,x3)
    axis equal
    legend('-DynamicLegend')
end

function [y] = integrator(ode,inicond,time,dt,p)

    y = NaN(length(time),length(inicond));
    y(1,:) = inicond;

  for j=2:length(time)
      k1 = dt*ode(y(j-1,:),p);
      k2 = dt*ode(y(j-1,:)+k1/2,p);
      k3 = dt*ode(y(j-1,:)+k2/2,p);
      k4 = dt*ode(y(j-1,:)+k3,p);
      y(j,:) = y(j-1,:) + k1/6+k2/3+k3/3+k4/6;
  end   
end



function [dydt] = ode(y,p)

dydt=NaN(1,6);

A = 3.25e-3;
B = 22e-3;
a = 100;
b = 50;
C = 135;
C1 = C;
C2 = 0.8*C;
C3 = 0.25*C;
C4 = 0.25*C;

dydt(1,1) = y(:,4);
dydt(1,4) = A*a*sigm(y(:,2)-y(:,3))-2*a*y(:,4)-(a^2)*y(:,1);
dydt(1,2) = y(:,5);
dydt(1,5) = A*a*(C2*sigm(C1*y(:,1)+p))-2*a*y(:,5)-(a^2)*y(:,2);
dydt(1,3) = y(:,6);
dydt(1,6) = B*b*C4*sigm(C3*y(:,1))-2*b*y(:,6)-(b^2)*y(:,3);
end

function S = sigm(u)
u0 = 6e-3;
e0 = 2.5;
r = 0.56e3;
S = 2*e0/(1+exp(r*(u0-u)));
end

Answer 1

There are a few comments you must observe to try to better understand the differences in the ode45 function and your non-built-in function integrator:

• The algorithm used in the ode45 function is NOT the classical 4th order Runge-Kutta method (1,2). It uses a different algorithm based on it. The ode45 help References section cite this paper (open archive) as the theory behind the implemented algorithm.

• Essentialy, this means the algorithm is an adaptive Runge-Kutta method. What the algorithm does is adjust the step size on each iteration. It determines the step size based on an estimate of the local truncation error.

• Knowing that, now you may understand better how the ode45 function works. If you pass as an argument the tspan with more than 2 elements, MATLAB will solve the differential equation with the steps chosen as by the algorithm and it will then interpolate the result in the given points:

If tspan has more than two elements [t0,t1,t2,...,tf], then the solver returns the solution evaluated at the given points. However, the solver does not step precisely to each point specified in tspan. Instead, the solver uses its own internal steps to compute the solution, then evaluates the solution at the requested points in tspan. The solutions produced at the specified points are of the same order of accuracy as the solutions computed at each internal step. Specifying several intermediate points has little effect on the efficiency of computation, but for large systems it can affect memory management.

• So, if your function has some "high degree of non regularity", ie, abruptly changes, discontinuities, high nonlinearities etc, for instance, the RK4 method works only if the time step is small enough.

• Bear in mind that all integrators can obtain the correct answer, as long as the step size is small enough. In practice, however, some functions are not well bevahed and the time step that must be used is so small that the computational time to solve the differential equation is prohibitively large .

Let's focus now on your particular problem. I will consider the case of p=-100 to illustrate the problem. I will solve it for a time span of 0 to 10, but your code solved it for 0 to 3 with the ode45 and 0 to 10 with the integrator.

Remember that, if you call the ode45 function without a left-hand argument, the ode opens a figure and plot the results by itself. If I don't specify a time step, notice the following:

p = -100;
time = [0 10];
initial = [0 0 0 0 0 0];
ode45(@ode,time,initial,[],p);

Observe the fast transient in the beginning. After about 0.1 sec, the solution is constant. Each o mark on the plot is one point utilized in the built-in ode45 function. The function determines this points.

We can call the ode45 function with a left-hand assignment and observe the time step determined by the ode:

[t,y] = ode45(@ode,time,initial,[],p);
plot(t(1:end-1),diff(t)); xlabel('t'); ylabel('\Deltat')

The time step utilized by the bult-in function is about 6.5e-3 in the permanent regime, but is smaller than that in the transient regime:

I = find(t<0.1);
plot(t(I(1:end-1)),diff(t(I))); xlabel('t'); ylabel('\Deltat')

The biggest time step in the transient regime is about 3e-3. It gets even smaller very close to the beginning (the abrupt changes, the discontinuities I mentioned earlier):

I = find(t<1e-4);
plot(t(I(1:end-1)),diff(t(I))); xlabel('t'); ylabel('\Deltat')

In summary, to solve this problem with the classical Runge-Kutta 4th order method you have to solve with a very small time step. This is not the recommended algorithm to solve your problem.