用欧拉法求解四维耦合系统

Question

我想用代码实现下面给出的系统，但是当我将它增加到 1500 次迭代时，我得到以下错误：

Warning (from warnings module):
  File "D:\python test files\sys1.py", line 16
    dy = c*x- x*z + w
RuntimeWarning: overflow encountered in double_scalars

Warning (from warnings module):
  File "D:\python test files\sys1.py", line 17
    dz = -b*z + x*y
RuntimeWarning: overflow encountered in double_scalars

Warning (from warnings module):
  File "D:\python test files\sys1.py", line 18
    du = -h*u - x*z
RuntimeWarning: overflow encountered in double_scalars

Warning (from warnings module):
  File "D:\python test files\sys1.py", line 42
    zs[i+1] = zs[i] + (dz * t)
RuntimeWarning: invalid value encountered in double_scalars

Warning (from warnings module):
  File "D:\python test files\sys1.py", line 15
    dx = a*(y-x) + u
RuntimeWarning: invalid value encountered in double_scalars

Warning (from warnings module):
  File "D:\python test files\sys1.py", line 19
    dw = k1*x - k2*y
RuntimeWarning: invalid value encountered in double_scalars

Warning (from warnings module):
  File "C:\Python27\lib\site-packages\mpl_toolkits\mplot3d\proj3d.py", line 156
    txs, tys, tzs = vecw[0]/w, vecw[1]/w, vecw[2]/w   
RuntimeWarning: invalid value encountered in divide

我的代码：

from __future__ import division
import numpy as np    
import math    
import random   
import matplotlib.pyplot as plt    
from mpl_toolkits.mplot3d import Axes3D        
# import pdb
# pdb.set_trace()

def sys1(x, y, z, u, w , a=10, b=8.0/3.0, c=28, k1=0.4, k2=8, h=-2):

    dx = a*(y-x) + u
    dy = c*x- x*z + w    
    dz = -b*z + x*y    
    du = -h*u - x*z    
    dw = k1*x - k2*y    
    return dx, dy, dz, du, dw

t = 0.01    
itera = 2500

# Need one more for the initial values

xs = np.empty((itera+1,))    
ys = np.empty((itera+1,))    
zs = np.empty((itera+1,))    
us = np.empty((itera+1,))    
ws = np.empty((itera+1,))

# Setting initial values

xs[0], ys[0], zs[0], us[0], ws[0] = (0.1, 0.1, 0.1, 0.1, 0.1)


# Stepping through "time".

for i in range(itera):    
# Derivatives of the X, Y, Z state
    dx, dy, dz, du, dw = sys1(xs[i], ys[i], zs[i], us[i], ws[i])     

    xs[i+1] = xs[i] + (dx * t)
    ys[i+1] = ys[i] + (dy * t)
    zs[i+1] = zs[i] + (dz * t)
    us[i+1] = us[i] + (du * t)
    ws[i+1] = ws[i] + (dw * t)

fig = plt.figure()

ax = fig.gca(projection='3d')
ax.plot(xs, ys, zs)
ax.set_xlabel("X Axis")
ax.set_ylabel("Y Axis")
ax.set_zlabel("Z Axis")
ax.set_title("Lorenz Attractor")
plt.show()

Answer 1

您正在尝试使用一个著名的简单数值方案来模拟一个著名的敏感非线性微分方程系统（实际上是a famously sensitive system 的增强 版本）。您的解决方案在给定的时间步长上出现分歧，这体现在您的解决方案值首先变为inf（您没有注意到），然后是nan（您仍然没有注意到），最后是@ 987654330@ 在你的无穷大附近玩转时产生除以零。

这是你的输出：

注意坐标轴的限制比例：都在 1e100 以上。这可能告诉你你有无穷无尽的东西。

好消息是，您拥有的同一个程序可以通过减少时间步长来提供合理的输出，这应该始终是您使用欧拉等一阶方法的第一个猜测（尤其是您从 MATLAB实现你的算法是正确的）。

使用t=0.001; itera=25000（左）和t=0.0001; itera=250000（右）的示例输出：

首先，请注意，这两个图是相当合理的，并且明显是有限的。其次，请注意，这两条轨迹虽然具有大致相似的形状，但非常不同。如果你使用更长的迭代（我的意思是更长的总迭代t*itera），差异会逐渐变得更加明显，最终两条轨迹将完全分道扬镳。

您应该确保您了解您正在使用一种非常不准确的方法来绘制一个非常敏感（准确地说：混沌）系统的轨迹。即使使用非常准确的方法，您最终也会积累一些错误，并且您会偏离初始值问题的实际解决方案。你所能期望的只是画出吸引子的大致形状，围绕它的轨迹将不可避免地呈锯齿形。但我相当肯定这就是你的目标。