python - 使用EULER'S方法求解4D耦合系统

我想用代码实现下面给出的系统，但是当我将其增加到1500次迭代时，就会出现以下错误：Warning (from warnings module): File "D:\python test files\sys1.py", line 16 dy = c*x- x*z + wRuntimeWarning: overflow encountered in double_scalarsWarning (from warnings module): File "D:\python test files\sys1.py", line 17 dz = -b*z + x*yRuntimeWarning: overflow encountered in double_scalarsWarning (from warnings module): File "D:\python test files\sys1.py", line 18 du = -h*u - x*zRuntimeWarning: overflow encountered in double_scalarsWarning (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_scalarsWarning (from warnings module): File "D:\python test files\sys1.py", line 15 dx = a*(y-x) + uRuntimeWarning: invalid value encountered in double_scalarsWarning (from warnings module): File "D:\python test files\sys1.py", line 19 dw = k1*x - k2*yRuntimeWarning: invalid value encountered in double_scalarsWarning (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]/wRuntimeWarning: invalid value encountered in divide我的代码：from __future__ import divisionimport numpy as npimport mathimport randomimport matplotlib.pyplot as pltfrom 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, dwt = 0.01itera = 2500# Need one more for the initial valuesxs = 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 valuesxs[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() 最佳答案 您正在尝试使用一个著名的简单数值方案来模拟一个著名的敏感的非线性微分方程组系统（实际上是a famously sensitive system的增强版本）。您的解决方案在给定的时间步长处发散，这首先在您的解决方案值中显示为inf（您没有注意到），然后变为nan（您仍然没有注意到），最后变为在处理您的无穷大时会产生零除。这是您的原样输出：注意轴的极限比例尺：所有都超过1e100。这本可以告诉您，您有无处不在。好消息是，您所拥有的同一程序仅需减少时间步就可以提供合理的输出，这始终是您使用Euler这样的一阶方法（尤其是您从MATLAB实现中确信，算法正确）。使用Axes3D.plot（左）和t=0.001; itera=25000（右）的示例输出： 首先，请注意，这两个图是相当合理的，并且公然有限。其次，请注意，这两个轨迹虽然形状大致相似，但差别很大。如果您使用更长的迭代次数（我的意思是更长的总t=0.0001; itera=250000），则差异将逐渐变得更加明显，最终两条轨迹将完全分开。您应该确保您了解使用的是非常不准确的方法来绘制非常敏感（准确地说是混沌）系统的轨迹。即使采用非常精确的方法，您最终也会累积一些错误，并且会从实际解决方案偏离您的初值问题。您所希望的就是绘制出吸引子的大致形状，围绕它的轨迹将不可避免地呈锯齿状。但我相当确定这是您的目标。关于python - 使用EULER'S方法求解4D耦合系统，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/36151796/