不带包装的Python中使用Runge Kutta 4th Order解决Lorentz模型 [英] Solving the Lorentz model using Runge Kutta 4th Order in Python without a package

问题描述

我希望在没有软件包帮助的情况下用Python解决Lorentz模型，并且我的代码似乎不符合我的期望.我不知道为什么我没有得到预期的结果和洛伦兹吸引者.我猜的主要问题与如何分别存储x，y和z的各种值有关.以下是我的带有3D解图的Lorentz模型的Runge-Kutta 45代码:

I wish to solve the Lorentz model in Python without the help of a package and my codes seems not to work to my expectation. I do not know why I am not getting the expected results and Lorentz attractor. The main problem I guess is related to how to store the various values for the solution of x,y and z respectively.Below are my codes for the Runge-Kutta 45 for the Lorentz model with 3D plot of solutions:

import numpy as np
import matplotlib.pyplot as plt
#from scipy.integrate import odeint

#a) Defining the Runge-Kutta45 method

def fx(x,y,z,t):
    dxdt=sigma*(y-z)
    return dxdt

def fy(x,y,z,t):
    dydt=x*(rho-z)-y
    return dydt

def fz(x,y,z,t):
    dzdt=x*y-beta*z
    return dzdt


def RungeKutta45(x,y,z,fx,fy,fz,t,h):
    k1x,k1y,k1z=h*fx(x,y,z,t),h*fy(x,y,z,t),h*fz(x,y,z,t)
    k2x,k2y,k2z=h*fx(x+k1x/2,y+k1y/2,z+k1z/2,t+h/2),h*fy(x+k1x/2,y+k1y/2,z+k1z/2,t+h/2),h*fz(x+k1x/2,y+k1y/2,z+k1z/2,t+h/2)
    k3x,k3y,k3z=h*fx(x+k2x/2,y+k2y/2,z+k2z/2,t+h/2),h*fy(x+k2x/2,y+k2y/2,z+k2z/2,t+h/2),h*fz(x+k2x/2,y+k2y/2,z+k2z/2,t+h/2)
    k4x,k4y,k4z=h*fx(x+k3x,y+k3y,z+k3z,t+h),h*fy(x+k3x,y+k3y,z+k3z,t+h),h*fz(x+k3x,y+k3y,z+k3z,t+h)
    return x+(k1x+2*k2x+2*k3x+k4x)/6,y+(k1y+2*k2y+2*k3y+k4y)/6,z+(k1z+2*k2z+2*k3z+k4z)/6



sigma=10.
beta=8./3.
rho=28.
tIn=0.
tFin=10.
h=0.05
totalSteps=int(np.floor((tFin-tIn)/h))

t=np.zeros(totalSteps)
x=np.zeros(totalSteps) 
y=np.zeros(totalSteps) 
z=np.zeros(totalSteps)

for i in range(1, totalSteps):
    x[i-1]=1. #Initial condition
    y[i-1]=1. #Initial condition
    z[i-1]=1. #Initial condition
    t[0]=0. #Starting value of t
    t[i]=t[i-1]+h
    x,y,z=RungeKutta45(x,y,z,fx,fy,fz,t[i-1],h)


#Plotting solution

from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm

fig=plt.figure()
ax=fig.gca(projection='3d')
ax.plot(x,y,z,'r',label='Lorentz 3D Solution')
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('z')
ax.legend()

推荐答案

我更改了集成步骤(顺便说一句，经典的四阶Runge-Kutta，而不是自适应RK54)，以广泛使用python核心概念的列表和列表操作以减少定义计算的位置数量.没有错误可以纠正，但是我认为算法本身更加集中.

I changed the integration step (btw., classical 4th order Runge-Kutta, not the adaptive RK54) to use the python core concept of lists and list operations extensively to reduce the number of places where the computation is defined. There were no errors there to correct, but I think the algorithm itself is more concentrated.

您的系统中出现错误，将其更改为快速发散的系统. Lorentz系统具有fx = sigma*(y-x)时，您拥有fx = sigma*(y-z).

You had an error in the system that changed it into a system that rapidly diverges. You had fx = sigma*(y-z) while the Lorentz system has fx = sigma*(y-x).

接下来，您的主循环有一些奇怪的任务.在每个循环中，您首先将先前的坐标设置为初始条件，然后用应用于完整数组的RK步骤替换完整数组.我完全替换掉了，找到正确解决方案的步伐不小.

Next your main loop has some strange assignments. In every loop you first set the previous coordinates to the initial conditions and then replace the full arrays with the RK step applied to the full arrays. I replaced that completely, there are no small steps to a correct solution.

import numpy as np
import matplotlib.pyplot as plt
#from scipy.integrate import odeint


def fx(x,y,z,t): return sigma*(y-x)
def fy(x,y,z,t): return x*(rho-z)-y
def fz(x,y,z,t): return x*y-beta*z

#a) Defining the classical Runge-Kutta 4th order method

def RungeKutta45(x,y,z,fx,fy,fz,t,h):
    k1x,k1y,k1z = ( h*f(x,y,z,t) for f in (fx,fy,fz) )
    xs, ys,zs,ts = ( r+0.5*kr for r,kr in zip((x,y,z,t),(k1x,k1y,k1z,h)) )
    k2x,k2y,k2z = ( h*f(xs,ys,zs,ts) for f in (fx,fy,fz) )
    xs, ys,zs,ts = ( r+0.5*kr for r,kr in zip((x,y,z,t),(k2x,k2y,k2z,h)) )
    k3x,k3y,k3z = ( h*f(xs,ys,zs,ts) for f in (fx,fy,fz) )
    xs, ys,zs,ts = ( r+kr for r,kr in zip((x,y,z,t),(k3x,k3y,k3z,h)) )
    k4x,k4y,k4z  =( h*f(xs,ys,zs,ts) for f in (fx,fy,fz) )
    return (r+(k1r+2*k2r+2*k3r+k4r)/6 for r,k1r,k2r,k3r,k4r in 
            zip((x,y,z),(k1x,k1y,k1z),(k2x,k2y,k2z),(k3x,k3y,k3z),(k4x,k4y,k4z)))



sigma=10.
beta=8./3.
rho=28.
tIn=0.
tFin=10.
h=0.01
totalSteps=int(np.floor((tFin-tIn)/h))

t = totalSteps * [0.0]
x = totalSteps * [0.0]
y = totalSteps * [0.0]
z = totalSteps * [0.0]

x[0],y[0],z[0],t[0] = 1., 1., 1., 0.  #Initial condition
for i in range(1, totalSteps):
    x[i],y[i],z[i] = RungeKutta45(x[i-1],y[i-1],z[i-1], fx,fy,fz, t[i-1], h)

使用tFin = 40和h=0.01我得到了图像

看起来像洛伦兹吸引子的典型图像.

looking like the typical image of the Lorentz attractor.

这篇关于不带包装的Python中使用Runge Kutta 4th Order解决Lorentz模型的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！