Вы можете переформулировать любое скалярное ODE (обыкновенное дифференциальное уравнение) порядка n в форме Коши в ODE порядка 1. Единственное, что вы «платите» в этой операции, это то, что переменные второго ODE будут векторами, а не скалярные функции.
Позвольте мне привести пример с ODE порядка 2. Предположим, ваш ODE: y '' = F (x, y, y '). Затем вы можете заменить его на [y, y ']' = [y ', F (x, y, y')], где производная вектора должна пониматься по компонентам.
Давайте вернем ваш код и вместо использования Рунге-Кутты порядка 4 в качестве приблизительного решения вашего ODE, мы применим простую схему Эйлера.
from pylab import*
import matplotlib.pyplot as plt
# we are approximating the solution of y' = f(x,y) for x in [x_0, x_1] satisfying the Cauchy condition y(x_0) = y0
def f(x, y0):
return y0
# here f defines the equation y' = y
def explicit_euler(x0, x1, y0, N,):
# The following formula relates h and N
h = (x1 - x0)/(N+1)
xd = list()
yd = list()
xd.append(x0)
yd.append(y0)
for i in range (1,N+1) :
# We use the explicite Euler scheme y_{i+1} = y_i + h * f(x_i, y_i)
y = yd[-1] + h * f(xd[-1], yd[-1])
# you can replace the above scheme by any other (R-K 4 for example !)
x = xd[-1] + h
yd.append(y)
xd.append(x)
return xd, yd
N = 250
x1 = 5
x0 = 0
y0 = 1
# the only function which satisfies y(0) = 1 and y'=y is y(x)=exp(x).
xd, yd =explicit_euler(x0, x1, y0, N)
plt.plot(xd,yd)
plt.show()
# this plot has the right shape !
Обратите внимание, что вы можете заменить схему Эйлера на R-K 4 , которая обладает лучшими свойствами стабильности и сходимости.
Теперь предположим, что вы хотите решить ODE второго порядка, скажем, например: y '' = -y с начальными условиями y (0) = 1 и y '(0) = 0. Затем вам нужно преобразовать Ваша скалярная функция y превращается в вектор размером 2, как описано выше и в комментариях в приведенном ниже коде.
from pylab import*
import matplotlib.pyplot as plt
import numpy as np
# we are approximating the solution of y'' = f(x,y,y') for x in [x_0, x_1] satisfying the Cauchy condition of order 2:
# y(x_0) = y0 and y'(x_0) = y1
def f(x, y_d_0, y_d_1):
return -y_d_0
# here f defines the equation y'' = -y
def explicit_euler(x0, x1, y0, y1, N,):
# The following formula relates h and N
h = (x1 - x0)/(N+1)
xd = list()
yd = list()
xd.append(x0)
# to allow group operations in R^2, we use the numpy library
yd.append(np.array([y0, y1]))
for i in range (1,N+1) :
# We use the explicite Euler scheme y_{i+1} = y_i + h * f(x_i, y_i)
# remember that now, yd is a list of vectors
# the equivalent order 1 equation is [y, y']' = [y', f(x,y,y')]
y = yd[-1] + h * np.array([yd[-1][1], f(xd[-1], yd[-1][0], yd[-1][1])]) # vector of dimension 2
print(y)
# you can replace the above scheme by any other (R-K 4 for example !)
x = xd[-1] + h # vector of dimension 1
yd.append(y)
xd.append(x)
return xd, yd
x0 = 0
x1 = 30
y0 = 1
y1 = 0
# the only function satisfying y(0) = 1, y'(0) = 0 and y'' = -y is y(x) = cos(x)
N = 5000
xd, yd =explicit_euler(x0, x1, y0, y1, N)
# I only want the first variable of yd
yd_1 = list(map(lambda y: y[0], yd))
plt.plot(xd,yd_1)
plt.show()
