@humoshaa
Астроном

Пытаюсь проинтегрировать уравнение методом Рунге-Кутта на питоне, сама задача поставлена правильно, проблема в программе?

Всем привет!
Уже которые день мучаюсь с данной мне задачей, аж пришлось задать вопрос в форуме. Я не большой эксперт по питону, всего лишь ученик, поэтому прошу простить если вопрос глуповат.
Мне нужно выполнить небольшую задачу с интегрированием уравнения, методом Рунге-Кутта, сам метод понятен. Но проблема возникает в самом коде, который прикреплю ниже:

#данные нам значения
x = [0.5]
y = [1.2]
y1 = [0.1]

#количество точек
x_last = 10.5

#шаг
delta_x = 0.1
n = int((x_last - x[0])/delta_x)

#производная
def dy(x, y, y1):
    return (y1*x)/(y**2+2.0)



def fi_0(x, y, y1):
    global delta_x
    return (delta_x/2)*dy(x, y, y1)

def fi_1(x, y, y1):
    return (delta_x/2) * dy(x + delta_x/2, y + y1*(delta_x/2), y1 + fi_0(x, y, y1))

def fi_2(x, y, y1):
    return (delta_x/2) * dy(x + delta_x/2, y + y1*(delta_x/2) + fi_0(x, y, y1)*(delta_x/2), y1 + fi_1(x, y, y1))

def fi_3(x, y, y1):
    return (delta_x/2) * dy(x + delta_x, y + delta_x*(y1 + fi_1(x, y, y1)), y1 + 2 * fi_2(x, y, y1))

def y_n(x, y, y1):
    return y + delta_x * (y1 + 1/3*(fi_0(x, y, y1) + fi_1(x, y, y1) + fi_2(x, y, y1)))
def y_n1(x, y, y1):
    return y1 + 1/3 * (fi_0(x, y, y1) + 2*fi_1(x, y, y1) + 2*fi_2(x, y, y1) + fi_3(x, y, y1))


for i in range(n):
    x0 = x[-1]; y0 = y[-1]; y10 = y1[-1]
    x.append(x0 + delta_x)
    y.append(y0 + y_n(x0, y0, y10))
    y1.append(y10 + y_n1(x0, y0, y10))

#вывод
for xi, yi, y1i in zip(x, y, y1):
    print(f'{xi:.2f}\t{yi:.8f}\t{y1i:.8f}')


Как я понял проблема в этой части кода:
for i in range(n):
    x0 = x[-1]; y0 = y[-1]; y10 = y1[-1]


Если изменю значения внутри квадратных скобок y0 и y10 на нули, то ответ похож на результат который должен быть, но все еще далеко от идеала. В нынешнем состоянии код показывает мне ответ в виде огромных чисел 15-17го порядка, а результат не должен превышать 2.2

Буду примного благодарен за любую помощь.
  • Вопрос задан
  • 176 просмотров
Пригласить эксперта
Ваш ответ на вопрос

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

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