Решение краевой задачи (уравнение диффузии-реакции) с помощью scipy solve_bvp - PullRequest
Новая
       101

Решение краевой задачи (уравнение диффузии-реакции) с помощью scipy solve_bvp

1 голос
/ 25 января 2020

Я изо всех сил пытаюсь решить следующую краевую задачу 2-го порядка: 

y'' + 2/x*y' + k**2.0*F(y) = 0

y(x=1)=1,  y'(x=0)=0

F(y) = -y or F(y) = -y*exp(AB*(1-y)/(1+B(1-y))

Я почему-то не могу правильно установить граничные условия. Я определил функцию для F(y)=y и граничных условий следующим образом: 

def fun(x, y, p):
 k = p[0]
 return np.vstack((y[1], -2.0/x*y[1] + k**2.0*y[0]))

def bc(ya, yb, p):
     return np.array([ya[0], yb[0],ya[1]])

y[0,:] = 1
y[1,0] = 0
sol = solve_bvp(fun, bc, x, y, p=[40])

Результаты, которые я должен получить, безусловно неверны, изменение начальных условий не улучшает ситуацию. Я думаю, что моя проблема в некоторой степени связана с граничным условием с нулевым градиентом при x = 0. Кто-нибудь знает, что я делаю неправильно?

РЕДАКТИРОВАТЬ: Здесь MWE, который должен дать постоянное значение 1 для k = 0,01. но для ak = 5 значение при x = 0 должно быть прибл. 0,06: 

def fun(x, y, p):
     k = p[0]
     return np.vstack((y[1], -2.0/x*y[1] + k**2.0*y[0]))

def bc(ya, yb, p):
     return np.array([ya[0], yb[0]-1.0,yb[1]])

x = np.linspace(1e-3, 1, 100)
y = np.zeros((2, x.size))
y[0,:] = 1

from scipy.integrate import solve_bvp
sol = solve_bvp(fun, bc, x, y, p=[1000])
x_plot = np.linspace(0, 1, 100)
y_plot = sol.sol(x_plot)[0]
plt.figure()
plt.plot(x_plot, y_plot)

1 Ответ

0 голосов
/ 26 января 2020

Рассмотрим случай F(y)=y. Тогда легко увидеть, что базисными решениями этого линейного ОДУ являются sin(k*x)/x и cos(kx/x). Аналогично для F(y)=-y получается sinh(k*x)/x и cosh(k*x)/x. Это означает, что большинство решений имеют особенность при x=0. С такой особенностью практически невозможно справиться со стандартными решениями ODE. Нужно помочь решателю в сингулярности, обычная процедура применяется снова на некотором расстоянии от сингулярности.

Что вы можете сделать, это проанализировать ситуацию в x=0 и немного отойти. Вы получаете через предел разностных отношений 

y''(0) + 2*y''(0) + k^2*F(y(0)) = 0

, который позволяет вычислить квадратичный c полином Тейлора. Таким образом, решите задачу на [a, 1] с помощью решателя ODE, используя продолжение к сингулярности y(x)=y(0)-k**2/6*F(y(0))*x**2 на [0,a].

Граничное условие на x=a проще всего установить sh с обработкой y0=y(0) в качестве параметра. Функции ODE и B C тогда равны 

def ode(x,y,y0): return [ y[1], -2*y[1]/x - k**2*F(y[0]) ]
def bc(ya,yb,y0): y2 = -k**2*F(y0)/6; return [ ya[0] - y0 - y2*a**2, ya[1] - 2*y2*a, yb[0]-1 ]

В случаях, обсуждаемых в вопросе, это дает 

a = 1e-2
def F(y): return -y

for k in [0.01, 5]:
    res = solve_bvp(ode, bc, [a,1], [[1,1], [0,0]], p=[1], tol=1e-5)
    print(f'k={k}: {res.message}, y0={res.p[0]}, theory: {k/np.sinh(k)}')
    if res.success:
        y0 = res.p[0]
        x = np.linspace(a,1,61);
        plt.plot(x, res.sol(x)[0])
        plt.plot([0], [y0],'o', res.x, res.y[0],'+', ms=4)
        plt.title(f'k={k}'); plt.grid(); plt.show()

с результатом

k = 0,01: алгоритм сходится с требуемой точностью., y0 = 0,9999933335277726, теория: 0,99999833335277757

plot of solution

k = 5: алгоритм сходится к требуемая точность., y0 = 0,06738256929427147, теория: 0,06738252915294544

plot of solution

...