Inconsistent solutions for numerical integration of non-linear differential equations?

I want to study the dynamics of the drive damped pendulum, given by the ode:

,

with constants and the rhs term the driving force.

I tried to solve it with the following code:

clear all
clc
close all

theta=0; dtheta=0; % Initial angle and angular velocity
tspan=[0 40]; % Integration time interval in seconds
y0=[theta;dtheta]; 
opts = odeset('MaxStep',0.001,"RelTol",1e-13,"AbsTol",1e-10);
[t,y]=ode45(@forced_damped_pendulum,tspan,y0,opts);

f1=figure(1); zoom=0.5;
f1.Position = [0 50 zoom*1900 zoom*900];
hold on
plot(t,y(:,1),'r','linewidth',1) % Plots the angle \theta (not dtheta)
xlabel('t','fontSize',14);
ylabel('theta','fontSize',14);
axis([tspan(1) tspan(2) -60 60]) % Eventually change the y-axis here
title('Chaotic Pendulum','fontsize',14)

function dydt=forced_damped_pendulum(t,y)
omega= 2*pi; omega0=(3/2)*omega; beta= omega0/4; gamma=1.16;
dydt=[y(2); gamma*omega0^2*cos(omega*t)-2*beta*y(2)-omega0^2*sin(y(1))];
end

The parameters have been chose to obtain chaotic dynamics, and thereby, the curve is higly irregular, as it should.

However, the problem is that whatever parameter options I use (MaxStep, RelTol, AbsTol, etc) and whatever ode (45, 78, 89, etc.) the result is always completely different. At best I get convergent results for few seconds (not much more than 20 second) but after that all solutions diverge form each other. Thus, the integration is not reliable. Not knowing what the "real" solution is, I can't see how to obtain a better numerical integration. Any thoughts, comments, reccomendetions?

[Torsten]

Didn't you answer your own question:

The parameters have been chose to obtain chaotic dynamics, and thereby, the curve is higly irregular, as it should.

?

Don't you think an integration method should reflect chaotic behaviour of the solution ?

[Marco]

Yes, but the same chaotic behavior, not a different one for whatever option and/or ode method.

[Torsten]

Yes, but the same chaotic behavior, not a different one for whatever option and/or ode method.

Chaotic means: a different one for whatever option and/or ode method.

There exists a unique solution, but we are not able to find it numerically.

Answer 1

[William Rose]

@Marco,

[Edit: Change "higher level solvers are able to achieve..." to "higher order solvers may achieve...", because the higher order solvers do not always use longer steps.]

Thanks for an interesting post! We are seeing a version of the butterfly effect in this chaotic system. The tiny differences between solvers cause the solutions to diverge. The higher order solvers are not necessarily more accurate, since stepsize is variable. The higher level solvers may achieve the specified accuracy requirtements with larger steps, i.e. with fewer function evaluations.

You said you want to know the real solution, so you can obtain better numerical integration. The solution with the highest level of accuracy which is computationally affordable to you (i.e. which can be done in a reasonable amount of time) is your best estimate of the real solution. If you tested these numerical solutions by comparison to experiment, I think you would find that experimental results for this system also diverge, due to the impossibility of doing exactly the same experiment twice.

Here is a plot of the results from your system, integrated with different solvers and accuracy options. Abbreviations "RT=-13, AT=-10, MS=-3" indicate relative tolerance=10^-13, absolute tolerance=10^-10, max step size=10^-3. If not specified, default options are used.

[Marco]

Thanks so far. Yes, these are the divergences that I described above. Of course, being a chaotic system, it is higly sensitive to intial conditions and to integration methods. But I hoped the latter could be minimized somehow. And that was precisley my intention: That of measuring how a real system is senstive to the initial conditions and to added noise. However, as I understand this, it is simply impossible to make a decently precise integration (say something converging to the "real" soilution for at least 100 seconds) becasue, not only we will never know what the exact solution (or at least almost exact in the 100 sec. interval) is, but no numerical solver will be ever able to find it anayway. Any attempt to measure such sensitivity makes no sense because we will never know how much the divergence is numerical and how much is due to the difference in inital conditions. Is that right? If so this is quite diaappointing, and I will have to give up everything. I suspect that there is a lot of nonsense in the scientific litarature about chaotic dynamics that presents similar integrations for much longer integration times.

[William Rose]

@Marco,

You write: "Any attempt to measure such sensitivity makes no sense because we will never know how much the divergence is numerical and how much is due to the difference in inital conditions. Is that right?"

I don't think that is right. You can use numerical integration to estimate the Lyapunov exponent, i.e. to estimate how much of the divergence is due to the difference in initial conditions. See here and here, for example. And in the example of using different solvers, explored above, the initial conditions are identical, so the divergence can be atributed to errors of numerical integration.

Others on this site know a lot more than me about chaotic dynamics.

Answer 2

[Marco]

But if every integration method gives different results, and different divergences, then they will also give different Lyapunov exponents.

[William Rose]

@Marco,

"if every integration method gives different results, and different divergences, then they will also give different Lyapunov exponents"

It is not obvious to me that that is true. Perhaps solutions associated with different initial conditions diverge faster than the divergence due to different solvers. If that is the case, then you may get consistent results for the estimate of the L.E., with different solvers.

If this were my interest, I would write code to implement the published methods for computing the Lyapounov exponent, or the Lyapunov spectrum, of a continuous system. I would check that my code gives the same result as the published papers, when analyzing those systems. Then I would apply that code to the nonlinear driven pendulum, with different solvers, to see if the computed LEs differ.