Константы Рунге-Кутты, расходящиеся для системы Лоренца?

Question

Я пытаюсь решить систему Лоренца, используя метод Рунге Кутты 4-го порядка, где

dx/dt=a*(y-x)
dy/dt=x(b-z)-y
dx/dt=x*y-c*z

Поскольку эта система не зависит от времени, возможно, эту часть итерации игнорировать, поэтому я просто имею dX=F(x,y,z)

def func(x0):
    a=10
    b=38.63
    c=8/3
    fx=a*(x0[1]-x0[0])
    fy=x0[0]*(b-x0[2])-x0[1]
    fz=x0[0]*x0[1]-c*x0[2]
    return np.array([fx,fy,fz])

def kcontants(f,h,x0):
    k0=h*f(x0)
    k1=h*f(f(x0)+k0/2)
    k2=h*f(f(x0)+k1/2)
    k3=h*f(f(x0)+k2)
    #note returned K is a matrix
    return np.array([k0,k1,k2,k3])

x0=np.array([-8,8,27])
h=0.001

t=np.arange(0,50,h)
result=np.zeros([len(t),3])

for time in range(len(t)):
    if time==0:
        k=kcontants(func,h,x0)
        result[time]=func(x0)+(1/6)*(k[0]+2*k[1]+2*k[2]+k[3])
    else:
        k=kcontants(func,h,result[time-1])
        result[time]=result[time-1]+(1/6)*(k[0]+2*k[1]+2*k[2]+k[3])

Результатом должны быть атракторы Лоренца, однако мой код расходится вокруг пятой итерации, и это потому, что константы, которые я создаю в kconstants сделайте, однако я проверил, и я почти уверен, что внедрение Rge Kutta не виновато... (по крайней мере, я думаю)

редактировать:

Нашел похожий пост, пока не могу понять, что я делаю не так

Answer 1

У вас есть дополнительный звонок f(x0) в расчете k1, k2 а также k3, Изменить функцию kcontants в

def kcontants(f,h,x0):
    k0=h*f(x0)
    k1=h*f(x0 + k0/2)
    k2=h*f(x0 + k1/2)
    k3=h*f(x0 + k2)
    #note returned K is a matrix
    return np.array([k0,k1,k2,k3])

Answer 2

Вы смотрели на разные начальные значения для вашего расчета? Имеют ли смысл те, которые вы выбрали? Т.е. они физические? Из прошлого опыта с rk вы можете иногда получить очень запутанные результаты, если выберете глупые стартовые параметры.

Answer 3

Спокойной ночи. Это и версия, которую я сделал с помощью интегратора scipy edo, scipy.integrate.odeint .

      # Author      : Carlos Eduardo da Silva Lima
# Theme       : Movement of a Plant around a fixed star
# Language    : Python
# date        : 11/19/2022
# Environment : Google Colab

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint
from scipy.optimize import root
from scipy.linalg import eig
from mpl_toolkits.mplot3d import Axes3D

##################################
# Condições inicial e parãmetros #
##################################
t_inicial = 0
t_final = 100
N = 10000
h = 1e-3

x_0 = 1.0
y_0 = 1.0
z_0 = 1.0

#####################
# Equação de Lorenz #
#####################
def Lorenz(r,t,sigma,rho,beta):

  x = r[0]; y = r[1]; z = r[2]

  edo1 = sigma*(y-x)
  edo2 = x*(rho-z)-y
  edo3 = x*y-beta*z
  return np.array([edo1,edo2,edo3])

t = np.linspace(t_inicial,t_final,N)
r_0 = np.array([x_0,y_0,z_0])
#sol = odeint(Lorenz,r_0,t,rtol=1e-6,args = (10,28,8/3))
sol = odeint(Lorenz, r_0, t, args=(10,28,8/3), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None, rtol=1e-9, atol=1e-9, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)

'''x = sol[:,0]
y = sol[:,1]
z = sol[:,2]'''
x, y, z = sol.T

# Plot
plt.style.use('dark_background')
ax = plt.figure(figsize = (10,10)).add_subplot(projection='3d')
ax.plot(x,y,z,'m-',lw=0.5, linewidth = 1.5)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Atrator de Lorenz")
plt.show()

Во второй части я моделирую две системы Лоренца, чтобы проверить чувствительные зависимости систем от начальных условий. Во второй системе я добавляю определенное количество eps =1e-3 к начальным условиям x(t0), y(t0) и z(t0).

      # Depedência com as condições iniciais
eps = 1e-3
r_0_eps = np.array([x_0+eps,y_0+eps,z_0+eps])
sol_eps = odeint(Lorenz, r_0_eps, t, args=(10,28,8/3), Dfun=None, col_deriv=0, full_output=0, ml=None, mu=None,
             rtol=1e-9, atol=1e-9, tcrit=None, h0=0.0, hmax=0.0, hmin=0.0, ixpr=0, mxstep=0, mxhnil=0, mxordn=12, mxords=5, printmessg=0, tfirst=False)

'''x_eps = sol_eps[:,0]
y_eps = sol_eps[:,1]
z_eps = sol_eps[:,2]'''
x_eps, y_eps, z_eps = sol_eps.T

# Plot
plt.style.use('dark_background')
ax = plt.figure(figsize = (10,10)).add_subplot(projection='3d')
ax.plot(x,y,z,'r-',lw=1.5)
ax.plot(x_eps,y_eps,z_eps,'b-.',lw=1.1)
ax.set_xlabel("X")
ax.set_ylabel("Y")
ax.set_zlabel("Z")
ax.set_title("Lorenz Attractor")
plt.show()

Надеюсь, я помог, увидимся :).