Как я могу остановить мой метод Runge-Kutta2 (Heun) от взрыва? - PullRequest
Новая
       96

Как я могу остановить мой метод Runge-Kutta2 (Heun) от взрыва?

1 голос
/ 23 февраля 2020

В настоящее время я пытаюсь написать некоторый код python для решения произвольной системы ODE первого порядка, используя общий явный метод Рунге-Кутты, определяемый значениями alpha, gamma (оба вектора измерения m) и beta (нижний tri angular матрица измерения mxm) таблицы Мясника, переданная пользователем. Мой код, кажется, работает для отдельных ODE, протестировав его на нескольких различных примерах, но я изо всех сил пытаюсь обобщить мой код для векторных ODE (то есть систем).

В частности, я пытаюсь решить ODE генератора Ван-дер-Поля (преобразованный в систему первого порядка), используя метод Хеуна, определенный значениями таблицы Батчера, приведенными в моем коде, но я получаю ошибки

  • "RuntimeWarning: переполнение, обнаруженное в double_scalars f = lambda t,u: np.array(... etc)" и
  • "RuntimeWarning: недопустимое значение, обнаруженное при добавлении kvec[i] = f(t+alpha[i]*h,y+h*sum)"

затем мой вектор решения, который явно взрывается. Обратите внимание, что закомментированный код, приведенный ниже, является одним из примеров отдельных ODE, которые я пробовал и решается правильно. Может ли кто-нибудь помочь, пожалуйста? Вот мой код: 

import numpy as np 

def rk(t,y,h,f,alpha,beta,gamma):
    '''Runga Kutta iteration'''
    return y + h*phi(t,y,h,f,alpha,beta,gamma)

def phi(t,y,h,f,alpha,beta,gamma):
    '''Phi function for the Runga Kutta iteration'''
    m = len(alpha)
    count = np.zeros(len(f(t,y)))
    kvec = k(t,y,h,f,alpha,beta,gamma)
    for i in range(1,m+1):
        count = count + gamma[i-1]*kvec[i-1]
    return count

def k(t,y,h,f,alpha,beta,gamma):
    '''returning a vector containing each step k_{i} in the m step Runga Kutta method'''
    m = len(alpha)
    kvec = np.zeros((m,len(f(t,y))))
    kvec[0] = f(t,y)
    for i in range(1,m):
        sum = np.zeros(len(f(t,y)))
        for l in range(1,i+1):
            sum = sum + beta[i][l-1]*kvec[l-1]
        kvec[i] = f(t+alpha[i]*h,y+h*sum)
    return kvec

def timeLoop(y0,N,f,alpha,beta,gamma,h,rk):
    '''function that loops through time using the RK method'''
    t  = np.zeros([N+1])
    y  = np.zeros([N+1,len(y0)])
    y[0] = y0
    t[0] = 0
    for i in range(1,N+1):
        y[i] = rk(t[i-1],y[i-1], h, f,alpha,beta,gamma)
        t[i] = t[i-1]+h
    return t,y

#################################################################

'''f  = lambda t,y: (c-y)**2
Y  = lambda t: np.array([(1+t*c*(c-1))/(1+t*(c-1))])
h0 = 1
c = 1.5
T = 10
alpha = np.array([0,1])
gamma = np.array([0.5,0.5])
beta = np.array([[0,0],[1,0]])
eff_rk   = compute(h0,Y(0),T,f,alpha,beta,gamma,rk, Y,11)'''

#constants
mu = 100
T = 1000
h = 0.01
N = int(T/h)

#initial conditions
y0 = 0.02
d0 = 0
init = np.array([y0,d0])

#Butcher Tableau for Heun's method
alpha = np.array([0,1])
gamma = np.array([0.5,0.5])
beta = np.array([[0,0],[1,0]])

#rhs of the ode system
f = lambda t,u: np.array([u[1],mu*(1-u[0]**2)*u[1]-u[0]])

#solving the system
time, sol = timeLoop(init,N,f,alpha,beta,gamma,h,rk)
print(sol)

1 Ответ

2 голосов
/ 23 февраля 2020

Ваш размер шага недостаточно мал. Генератор Ван-дер-Поля с mu=100 - это быстродействующая медленная система с очень резкими поворотами при переключении режимов, поэтому довольно жесткая. При использовании явных методов это требует малых размеров шагов, наименьший разумный размер шагов составляет от 1e-5 до 1e-6. Вы получаете решение по предельному циклу уже для h=0.001, с результирующими скоростями до 150.

Вы можете уменьшить часть этой жесткости, используя другую переменную скорость / импульс. В уравнении 

x'' - mu*(1-x^2)*x' + x = 0

вы можете объединить первые два члена в производную, 

mu*v = x' - mu*(1-x^2/3)*x

, чтобы 

x' = mu*(v+(1-x^2/3)*x)
v' = -x/mu

Второе уравнение теперь было равномерно медленным близко к предельному циклу, в то время как первый имеет длинные относительно прямые прыжки, когда v покидает куби c v=x^3/3-x.

Это прекрасно интегрируется с оригинальным h=0.01, сохраняя решение внутри коробки [-3,3]x[-2,2].

...