洛伦兹吸引子与Runge-Kutta蟒蛇 [英] Lorenz attractor with Runge-Kutta python

问题描述

1. 您好，我必须编写一个python函数以使用Runge-Kutta 2cond级来求解Lorenz微分方程

sigma = 10，r = 28，b = 8/3

初始条件为(x，y，z)=(0,1,0)

这是我编写的代码，但是它抛出一个错误，提示在 double_scalars 中遇到溢出，我看不出程序有什么问题

从pylab导入

  *def runge_4(r0，a，b，n，f1，f2，f3):def f(r，t):x = r [0]y = r [1]z = r [2]fx = f1(x，y，z，t)fy = f2(x，y，z，t)fz = f3(x，y，z，t)返回数组([fx，fy，fz]，float)h =(b-a)/nlista_t =范围(a，b，h)打印(lista_t)X，Y，Z = []，[]，[]对于lista_t中的t:k1 = h * f(r0，t)print("k1 ="，k1)k2 = h * f(r0 + 0.5 * k1，t + 0.5 * h)print("k2 ="，k2)k3 = h * f(r0 + 0.5 * k2，t + 0.5 * h)print("k3 ="，k3)k4 = h * f(r0 + k3，t + h)print("k4 ="，k4)r0 + =(k1 + 2 * k2 + 2 * k3 + k4)/float(6)打印(r0)X.append(r0 [0])Y.append(r0 [1])Z.append(r0 [2])返回数组([X，Y，Z])def f1(x，y，z，t):返回10 *(y-x)def f2(x，y，z，t):返回28 * x-y-x * zdef f3(x，y，z，t):返回x * y-(8.0/3.0)* z#然后我运行r0 = [1,1,1]runge_4(r0,1,50,20，f1，f2，f3) 

解决方案

以数字方式求解微分方程可能具有挑战性.如果您选择太大的步长，则解决方案将积累很多错误，甚至可能变得不稳定，就像您的情况一样.

您应该大大减小步长( h )，或者只使用 scipy 提供的自适应Runge Kutta方法:

来自numpy导入数组的

 从scipy.integrate导入solve_ivp导入pylab从mpl_toolkits导入mplot3ddef func(t，r):x，y，z = rfx = 10 *(y-x)fy = 28 * x-y-x * zfz = x * y-(8.0/3.0)* z返回数组([fx，fy，fz]，float)#我运行它r0 = [0，1，0]溶胶= resolve_ivp(func，[0，50]，r0，t_eval = linspace(0，50，5000))#并绘制无花果= pylab.figure()斧= pylab.axes(projection ="3d")ax.plot3D(sol.y [0 ,:]，sol.y [1 ,:]，sol.y [2 ,:]，'blue')pylab.show() 

此求解器使用4阶和5阶Runge Kutta组合，并通过调整步长来控制它们之间的偏差.在此处查看更多用法文档:

1. Hello I have to program a python function to solve Lorenz differential equations using Runge-Kutta 2cond grade

sigma=10, r=28 and b=8/3

with initial conditions (x,y,z)=(0,1,0)

this is the code i wrote, but it throws me an error saying overflow encountered in double_scalars, and I don't see what is wrong with the program

from pylab import *
def runge_4(r0,a,b,n,f1,f2,f3):
    def f(r,t):
        x=r[0]
        y=r[1]
        z=r[2]
        fx=f1(x,y,z,t)
        fy=f2(x,y,z,t)
        fz=f3(x,y,z,t)
        return array([fx,fy,fz],float)
    h=(b-a)/n
    lista_t=arange(a,b,h)
    print(lista_t)
    X,Y,Z=[],[],[]
    for t in lista_t:
        k1=h*f(r0,t)
        print("k1=",k1)
        k2=h*f(r0+0.5*k1,t+0.5*h)
        print("k2=",k2)
        k3=h*f(r0+0.5*k2,t+0.5*h)
        print("k3=",k3)
        k4=h*f(r0+k3,t+h)
        print("k4=",k4)
        r0+=(k1+2*k2+2*k3+k4)/float(6)
        print(r0)
        X.append(r0[0])
        Y.append(r0[1])
        Z.append(r0[2])
    return array([X,Y,Z])

def f1(x,y,z,t):
    return 10*(y-x)
def f2(x,y,z,t):
    return 28*x-y-x*z
def f3(x,y,z,t):
    return x*y-(8.0/3.0)*z
#and I run it
r0=[1,1,1]

runge_4(r0,1,50,20,f1,f2,f3)

解决方案

Solving differential equations numerically can be challenging. If you choose too high step sizes, the solution will accumulate high errors and can even become unstable, as in your case.

Either you should drastically reduce the step size (h) or just use the adaptive Runge Kutta method provided by scipy:

from numpy import array, linspace
from scipy.integrate import solve_ivp
import pylab
from mpl_toolkits import mplot3d

def func(t, r):
    x, y, z = r 
    fx = 10 * (y - x)
    fy = 28 * x - y - x * z
    fz = x * y - (8.0 / 3.0) * z
    return array([fx, fy, fz], float)

# and I run it
r0 = [0, 1, 0]
sol = solve_ivp(func, [0, 50], r0, t_eval=linspace(0, 50, 5000))

# and plot it
fig = pylab.figure()
ax = pylab.axes(projection="3d")
ax.plot3D(sol.y[0,:], sol.y[1,:], sol.y[2,:], 'blue')
pylab.show()

This solver uses 4th and 5th order Runge Kutta combination and controls the deviation between them by adapting the step size. See more usage documentation here: https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_ivp.html

这篇关于洛伦兹吸引子与Runge-Kutta蟒蛇的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！