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

4 Ansichten (letzte 30 Tage)
Marco am 17 Jan. 2024
Kommentiert: William Rose am 18 Jan. 2024
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?
##### 3 Kommentare1 älteren Kommentar anzeigen1 älteren Kommentar ausblenden
Marco am 17 Jan. 2024
Yes, but the same chaotic behavior, not a different one for whatever option and/or ode method.
Torsten am 17 Jan. 2024
Bearbeitet: Torsten am 17 Jan. 2024
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.

Melden Sie sich an, um zu kommentieren.

### Antworten (2)

William Rose am 17 Jan. 2024
Bearbeitet: William Rose am 17 Jan. 2024
[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.
##### 2 KommentareKeine anzeigenKeine ausblenden
Marco am 17 Jan. 2024
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 am 18 Jan. 2024
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.

Melden Sie sich an, um zu kommentieren.

Marco am 18 Jan. 2024
But if every integration method gives different results, and different divergences, then they will also give different Lyapunov exponents.
##### 1 Kommentar-1 ältere Kommentare anzeigen-1 ältere Kommentare ausblenden
William Rose am 18 Jan. 2024
"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.

Melden Sie sich an, um zu kommentieren.

### Kategorien

Mehr zu Matrix Computations finden Sie in Help Center und File Exchange

R2023b

### Community Treasure Hunt

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

Start Hunting!

Translated by