Inconsistent solutions for numerical integration of non-linear differential equations?
7 次查看(过去 30 天)
显示 更早的评论
I want to study the dynamics of the drive damped pendulum, given by the ode:
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1593071/image.png)
with
constants and the rhs term the driving force.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1593076/image.png)
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 个评论
回答(2 个)
William Rose
2024-1-17
编辑:William Rose
2024-1-17
[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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1593186/image.jpeg)
2 个评论
William Rose
2024-1-18
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.
Marco
2024-1-18
1 个评论
William Rose
2024-1-18
"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.
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Matrix Computations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!