Решение модели Лоренца с использованием Runge Kutta 4-го порядка в Python без пакета

Question

Я хочу решить модель Лоренца в Python без помощи пакета, и мои коды, похоже, не работают, как я ожидал. Я не знаю, почему я не получаю ожидаемых результатов и аттрактор Лоренца. Основная проблема, которую я предполагаю, связана с тем, как хранить различные значения для решения x,y и z соответственно. Ниже приведены мои коды для Runge-Kutta 45 для модели Лоренца с трехмерным графиком решений:

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()

Answer 1

Я изменил шаг интеграции (между прочим, классический Runge-Kutta 4-го порядка, а не адаптивный RK54), чтобы использовать основную концепцию Python для списков и операций со списками, чтобы уменьшить количество мест, где определяются вычисления. Там не было ошибок, чтобы исправить, но я думаю, что сам алгоритм более концентрированный.

У вас была ошибка в системе, которая превратила ее в систему, которая быстро расходится. Ты имел fx = sigma*(y-z) в то время как система Лоренца имеет fx = sigma*(y-x),

Далее ваш главный цикл имеет несколько странных назначений. В каждом цикле вы сначала устанавливаете предыдущие координаты в начальные условия, а затем заменяете полные массивы шагом RK, примененным к полным массивам. Я заменил это полностью, нет никаких маленьких шагов к правильному решению.

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 Я получаю изображение

выглядит как типичное изображение аттрактора Лоренца.