使用 EULER 方法求解 4D 耦合系统

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（您仍然没有注意到），最后Axes3D.plot 中的缩放在处理无穷大时产生除以零。

这是您的原样输出：

注意坐标轴的极限比例：都在1e100以上。这本可以告诉您周围有无穷大运行。

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

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

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

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

使用 EULER 方法求解 4D 耦合系统

Solving 4D coupled system by using EULER'S Method

python

python-2.7

numpy

matplotlib

ode