Как изменить реализацию этого метода Рунге-Кутты с решателя ода первого порядка на систему решателя ОДУ? - PullRequest
Новая
       94

Как изменить реализацию этого метода Рунге-Кутты с решателя ода первого порядка на систему решателя ОДУ?

0 голосов
/ 24 апреля 2019

Я реализовал 4-х шаговый метод Рунге-Кутты (k1..k4) ODE-решатель для функции u ′ (x) = f (x, u (x)) с начальным условием u (x0) = u0

Но он решает только ОДУ первого порядка.Как я могу изменить код, чтобы он также принимал ODE более высокого уровня? u ′ (x) = f (x, u (x)) с начальным условием u (x0) = u0 ?Все как векторы здесь не скаляры.

Вот код для моего решения, в котором используются ODE первого порядка: 

def euler(f, x, y, h):
    yn = y + h*f(x,y)
    xn = x + h
    return xn, yn

def rk4(f, x, y, h):
    k1 = h*f(x,y)
    k2 = h*f(x+(1/2)*h,y+(1/2)*k1) 
    k3 = h*f(x+(1/2)*h,y+(1/2)*k2)
    k4 = h*f(x+h,y+k3)
    yn = y + ((1/6)*(k1+(2*k2)+(2*k3)+k4))
    xn = x + h
    return xn, yn

def integrate(method, f, t0, y0, tend, h): # Depending on the 'method' (euler or rk4) this method solves the ode f)
    T = [t0]
    Y = [y0]
    t = t0
    y = y0 
    while (t < tend):
        h = min(h, tend-t)
        t, y = method(f, t, y, h)
        T.append(t)
        Y.append(y)
    return np.array(T), np.array(Y)

def f1(t, y):   # this is an 1st order ODE
    return (t*y) + t

# Usage example
xv, yv = integrate(method=rk4, f=f1, t0=0, y0=2, tend=1, h=0.5)
yv[-1] # Output: 6.15203857421875

Есть ли решение, позволяющее сделать эту работу также для данной системы ODE?(динамически, без знания входного ODE, если он 1-го или более высокого порядка), например, этот: u ′ ′ (x) + u (x) + u (x) ^ 3 = 0 с u (0) = u ′ (0) = 0 или этот: u′1 (x) = 98u1 (x) + 198u2 (x) и u′2 (x) = −99u1 (x) - 199u2 (x) с u1 (0) = 1 и u2(0) = 0

...