MisticX
@MisticX
Кратно обо мне

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

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

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

Дело в том, что уже на 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)


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

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

https://www.cyberforum.ru/cpp-beginners/thread1874... - отсюда я взял коэффициенты(последнее сообщение с кодом)
  • Вопрос задан
  • 100 просмотров
Решения вопроса 1
hint000
@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
Ответ написан
Пригласить эксперта
Ваш ответ на вопрос

Войдите, чтобы написать ответ

Войти через центр авторизации
Похожие вопросы