python - 使用 Lorenz 方程的运行时警告
问题描述
我正在尝试在 python 中做 Lorenz 方程(我正在关注 Mark Newman - Computational Physics(2012,CreateSpace Independent Publishing Platform)的练习 8.3)我已经得到了图形,一切看起来“正确”。这可能是一个数学问题,而不是一个真正的编程问题,但我在这里发布以确保。首先,这是我的代码:
from numpy import arange,array
import pylab
def f(v,t):
s=10
r=28
b=8/3
x= v[0]
y= v[1]
z= v[2]
fx= s*(y - x)
fy= r*x - y - x*z
fz= x*y - b*z
return array([fx,fy,fz],float)
def d(N):
a=0.0
b=50.0
h=(b-a)/N
r=array([0.0,1.0,0.0],float)
tpoints=arange(a,b,h)
xpoints= []
ypoints= []
zpoints= []
for t in tpoints:
xpoints.append(r[0])
ypoints.append(r[1])
zpoints.append(r[2])
k1 = h*f(r,t)
k2 = h*f(r+0.5*k1,t+0.5*h)
k3 = h*f(r+0.5*k2,t+0.5*h)
k4 = h*f(r+k3,t+h)
r += (k1+2*k2+2*k3+k4)*(1/6)
return tpoints,xpoints,ypoints,zpoints
for i in range (1,6):
N=10**i
pylab.plot(d(N)[0],d(N)[1], label=N)
pylab.xlabel("t")
pylab.ylabel("x(t)")
pylab.title("Gráficos x em função de t")
pylab.legend()
pylab.show()
pylab.plot(d(N)[0],d(N)[2], label=N)
pylab.xlabel("t")
pylab.ylabel("y(t)")
pylab.title("Gráficos y em função de t")
pylab.legend()
pylab.show()
pylab.plot(d(N)[0],d(N)[3], label=N)
pylab.xlabel("t")
pylab.ylabel("z(t)")
pylab.title("Gráficos z em função de t")
pylab.legend()
pylab.show()
pylab.plot(d(N)[1],d(N)[3], label=N)
pylab.xlabel("x")
pylab.ylabel("z(x)")
pylab.title("Gráficos z em função de x")
pylab.legend()
pylab.show()
这给了我解决问题的图表,我认为这是正确的。
当我去它i=1
给i=3
我这个错误:
我认为这与数学问题有关,但是当我搜索错误时,它给了我一些带有数组的东西。所以我正在检查它。
对于i
等于或大于 4,它没有给我带来任何问题。
解决方案
带有 RK4 的 Lorenz 系统需要大约或更少的步长0.05
,也就是说,N=10**4
对于您的构造,或更大。使用 Runge-Kutta python查看近似重复的Lorenz 吸引子。
对于较大的步长,即出现错误的情况,该方法将返回混沌结果,这些结果大多与系统的精确解以及与之相关的任何界限无关。因此,由于二次超线性项,浮点溢出的发散是可能的。
推荐阅读
- spring - Static resource mapping in Spring3
- javascript - Javascript 滚动和验证表单
- cloud - 如何在juju的ocata openstack中部署多区域?
- python - 同比 年初至今 百分比变化
- java - 使用 COUNT(*) 从 netbean 中的多个表中计算多行
- xamarin.forms - Xamarin Forms UWP 使用参数启动应用程序
- git - 在提交历史中附加一个 git 提交
- angular - Angular 两个按钮提交表单但用途不同
- php - Telegram Bot php 键盘数组按钮动态
- java - 构建抽象类的实例是什么意思?