Как исправить ошибку при вычислениях методом Рунге-Кутта (4 порядка)?

Question

[Максим Ко]

Доброго времени суток!

Появилась проблема с вычислением значений для построения графика аттрактора Лоренца.

Дело в том, что уже на 2 - 3 итерации начинают появляться очень большие значения коэффициентов.

Код:

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

sigma = 10
r = 28
b = 8 / 3

def LorenX(x,y):
    return sigma * (x - y)
def LorenY(x,y,z):
    return x * (r - z) - y
def LorenZ(x,y,z):
    return x * y - b * z

def rg4(x0,y0,z0,dt = 0.0001):
    x = x0
    y = y0
    z = z0
    arrx = [0]
    arry = [0]
    arrz = [0]

    fig = plt.figure()
    ax = fig.add_subplot(111, projection='3d')
    for i in range(10000):

        arrx.append(x)
        arry.append(y)
        arrz.append(z)

        k1 = LorenX(x,y)
        m1 = LorenY(x,y,z)
        p1 = LorenZ(x,y,z)

        k2 = LorenX(x + k1/2,y + m1/2)
        m2 = LorenY(x + k1/2,y + m1/2,z + p1/2)
        p2 = LorenZ(x + k1/2,y + m1/2,z + p1/2)

        k3 = LorenX(x + k2/2,y + m2/2)
        m3 = LorenY(x + k2/2,y + m2/2,z + p2/2)
        p3 = LorenZ(x + k2/2,y + m2/2,z + p2/2)

        k4 = LorenX(x + k3,y + m3)
        m4 = LorenY(x + k3,y + m3,z + p3)
        p4 = LorenZ(x + k3,y + m3,z + p3)

        x += (dt * (k1 + 2*k2 + 2*k3 + k4) / 6)
        y += (dt * (m1 + 2*m2 + 2*m3 + m4) / 6)
        z += (dt * (p1 + 2*p2 + 2*p3 + p4) / 6)
    ax.plot(arrx, arry, arrz)
    plt.show()

rg4(1,1,1)

Здесь приведены первые две итерации (значения x,y,z соответственно)

Pyplot строит такую линию и три оси с неадекватными значениями:


https://www.cyberforum.ru/cpp-beginners/thread1874... - отсюда я взял коэффициенты(последнее сообщение с кодом)

Answer 1

[hint000]

Во-первых:
# b = 8 / 3
# получилось целочисленное деление, b=2
b = 8.0/3
Во-вторых:
def LorenX(x,y):
# return sigma * (x - y)
return sigma * (y - x)
В-третьих:

...отсюда я взял коэффициенты(последнее сообщение с кодом)

Неправильная копипаста же. Что-то странное с dt вы сделали. Я оттуда же скопипастил без отсебятины:

#
        k1 = dt*LorenX(x,y)
        m1 = dt*LorenY(x,y,z)
        p1 = dt*LorenZ(x,y,z)

        k2 = dt*LorenX(x + k1/2,y + m1/2)
        m2 = dt*LorenY(x + k1/2,y + m1/2,z + p1/2)
        p2 = dt*LorenZ(x + k1/2,y + m1/2,z + p1/2)

        k3 = dt*LorenX(x + k1/2,y + m2/2)
        m3 = dt*LorenY(x + k1/2,y + m2/2,z + p2/2)
        p3 = dt*LorenZ(x + k1/2,y + m2/2,z + p2/2)

        k4 = dt*LorenX(x + k1,y + m3)
        m4 = dt*LorenY(x + k1,y + m3,z + p3)
        p4 = dt*LorenZ(x + k1,y + m3,z + p3)

        x += (k1 + 2*k2 + 2*k3 + k4) / 6
        y += (m1 + 2*m2 + 2*m3 + m4) / 6
        z += (p1 + 2*p2 + 2*p3 + p4) / 6

Comments

[Максим Ко]

Спасибо! Что-то я даже не обратил внимание коэффициенты и просто подумал, что вынес за скобки)

[Andy_U]

Во-первых:
# b = 8 / 3
# получилось целочисленное деление, b=2

Вы до сих пор используете второй питон?

D:\Temp>py
Python 3.10.2 (tags/v3.10.2:a58ebcc, Jan 17 2022, 14:12:15) [MSC v.1929 64 bit (AMD64)] on win32
Type "help", "copyright", "credits" or "license" for more information.
>>> b=8/3
>>> b
2.6666666666666665

[hint000]

Andy_U, я вообще не использую Питон и не знаю его. :) Просто в Линуксе и так "из коробки" есть и 2-ой (для совместимости, видимо) и 3-й.

username@host:~$ py
Command 'py' not found, but can be installed with:
sudo apt install pythonpy

username@host:~$ python
Python 2.7.18 (default, Feb 22 2022, 11:45:08) 
[GCC 11.2.0] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> 

username@host:~$ python3
Python 3.10.2 (main, Feb 26 2022, 08:21:15) [GCC 11.2.0] on linux
Type "help", "copyright", "credits" or "license" for more information.
>>>

[Andy_U]

hint000, Что не используете, извините, заметно. Писать свою функцию для решения обычного дифф.уравнения, вместо использования scipy это несколько странно.

[hint000]

Andy_U, вы меня перепутали с автором вопроса. :) Я никаких своих функций не писал, я лишь исправил ошибки.

[Andy_U]

Не перепутал. Никому бы из плотно работающих с питоном не пришло бы в голову использовать свой велосипед, где даже numpy не используется (код стал бы в 3 раза короче). Где автоматический выбор шага?