У меня есть два алгоритма для задачи численного дифференциального уравнения, один называется методом Эйлера, а второй называется Рунге Кутта второго порядка (RK2). По существу метод Эйлера и RK2 приближают решение дифференциальных уравнений. Единственное отличие состоит в том, что они используют разные формулы (Эйлер использует производную первого порядка ряда Тейлора, в то время как RK2 является производной второго порядка ряда Тейлора).
Я пытаюсь исправить некоторый написанный мной код чтобы он мог вернуть следующее решение, 
Однако мой код не возвращает правильные значения, когда дело доходит до RK2, но возвращает следующее,
Обратите внимание, что в моем решении размер шага h обозначается как dt. Я предоставил код, который использовал для создания этого ниже, а затем приведен числовой c пример метода Рунге Кутты второго порядка, который работает численно. Мне интересно показать, что с RK2 сходимость происходит быстрее, чем по методу Эйлера.
import matplotlib.pyplot as plt
import numpy as np
from math import exp # exponential function
dy = lambda x, y: x * y
f = lambda x: exp(x ** 2 / 2) # analytical solution function
x_final = 2
# analytical solution
x_a = np.arange(0, x_final, 0.01)
y_a = np.zeros(len(x_a))
for i in range(len(x_a)):
y_a[i] = f(x_a[i])
plt.plot(x_a, y_a, label="analytical")
# Container for step sizes dt /dt
dt = 0.5
x = 0
y = 1
print("dt = " + str(dt))
print("x \t\t y (Euler) \t y (analytical)")
print("%f \t %f \t %f" % (x, y, f(x)))
n = int((x_final - x) / dt)
x_n = np.zeros(n + 1)
y_n = np.zeros(n + 1)
x_n[0] = x
y_n[0] = y
#Plot Euler's method
for i in range(n):
y += dy(x, y) * dt
x += dt
print("%f \t %f \t %f" % (x, y, f(x)))
x_n[i + 1] = x
y_n[i + 1] = y
plt.plot(x_n, y_n, "x-", label="Euler dt=" + str(dt))
###################33
# Runge-Kutta's method 2'nd order (RK2)
x = 0
y = 1
print("dt = " + str(dt))
print("x \t\t y (rk2) \t y (analytical)")
print("%f \t %f \t %f" % (x, y, f(x)))
n = int((x_final - x) / dt)
x_n = np.zeros(n + 1)
y_n = np.zeros(n + 1)
x_n[0] = x
y_n[0] = y
# Plot the RK2
for z in range(n):
K1 = dt*dy(x,y) # Step 1
K2 = dt*dy(x+dt/2,y+K1/K2) # Step 2
y += K2 # Step 3
x += dt
print("%f \t %f \t %f" % (x, y, f(x)))
x_n[i + 1] = x
y_n[i + 1] = y
plt.plot(x_n, y_n, "x-", label="RK2 dt=" + str(dt))
plt.title("Euler's method ")
plt.xlabel("x")
plt.ylabel("y")
plt.legend()
plt.show()
Это код с цифрой c, который работает для RK2. Я поместил шаг 1,2 и 3 в изменения метода Эйлера, который позволяет нам создать RK2.
from math import exp
dy = lambda x,y: x*y
f = lambda x: exp(x**2/2)
x = 0
xn = 2
y = 1
dt = 0.5
n = int((xn)/dt)
print ('x \t\t y (RK2) \t y (analytical)')
print ('%f \t %f \t %f'% (x,y,f(x)))
# main loop
for i in range(n):
K1 = dt*dy(x, y) # step 1
K2 = dt*dy(x + dt/2, y + K1/2) # step 2
y += K2 # step 3
x += dt
print ('%f \t %f \t %f'% (x,y,f(x)))
Извините, если я задал этот вопрос плохо. Я новичок в python, поэтому я все еще пытаюсь понять, как go решать эти типы проблем. Вкратце мой вопрос состоит в том, как я могу исправить вышеуказанную функцию, чтобы вычислить правильную оценку и затем отобразить ее на графике, используя matplotlib из python.