SciPy: solve_bvp Задача 2-й порядок Дифф. уравнение - PullRequest
Новая
       69

SciPy: solve_bvp Задача 2-й порядок Дифф. уравнение

0 голосов
/ 19 февраля 2020

Пытается решить разницу 2-го порядка. э. с двумя граничными условиями и ничего, что я пытаюсь сделать, кажется, не работает, и я не могу найти учебник, который включает все / подобные термины, которые есть в моем выражении, и, по крайней мере для меня, документация scipy не объясняет, как использовать Ясно решить. bvp.

У меня есть уравнение: y '' + 2 / r * y '= A * y + b * y ^ 3, где y - функция от r.

I' мы переписали его следующим образом:

y1 = y (r)

y2 = y1 '

, так что y2' = -2 / r * y2 + y1 (A + b * y1 ^ 2)

При условии, что y '(0) = 0, y (r = 10) = постоянная

A, b и постоянная известны.

У меня есть следующий код, но он, кажется, не работает, и, как уже было сказано, я немного смущен документацией, поэтому любая помощь будет оценена! 

def fun(x, y):
    return np.vstack((y[2], -2/x*y[1]+y[0]*(A+b*y[0]*y[0])))

def bc(ya, yb):
    return np.array([ya[2], yb[1]-constant])


x = np.linspace(1, 10, 10)
ya = np.zeros((3, x.size))
yb = np.zeros((3, x.size))

sol_1 = solve_bvp(fun, bc, x, ya)
sol_2 = solve_bvp(fun, bc, x, yb)

Спасибо!

========== РЕДАКТИРОВАТЬ ========================= Есть аналитическое решение, для которого я рассчитал Я просто смотрю, смогу ли я найти то же решение численно, я думаю, что главная проблема заключается в том, что Для решения есть две отдельные области: одна, где r R (out). Это приводит к двум различным решениям (действительным только в их соответствующей области) с условиями, что y_in (R) == y_out (R) и y_in '(R) == y_out' (R). Полное решение, состоящее из двух частей, где Радиус = 1, a = 99, b = 1 и постоянная = 1, y (inf) = постоянная

Из решения Лутца Лемана он принимает правильную форму (по крайней мере, для внутренней области, хотя и не в правильных масштабах).

Я просто не уверен, как вы будете go кодировать все уравнивающие решения, и я полагаю, даже получая их решения в первую очередь, хотя ответ Лутца - удивительный момент в правильном направлении. Спасибо

1 Ответ

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

Выпуски

  • Порядок уравнения равен 2, поэтому размерность вектора состояния равна 2, значение всегда y[0], производная y[1], есть нет y[2], вероятно, остаток от перевода из Matlab.

  • Также в граничных условиях нет ya[2], производное значение равно ya[1], а функция значение во втором - yb[0].

  • Исходное предположение должно иметь одинаковое число 2 компонентов состояния.

  • Почему Вы рассчитываете два решения с одинаковыми данными?

  • Примечание: нет необходимости преобразовывать возвращаемые значения в типы numpy, решатель проверяет и конвертирует в любом случае.

Каркас BVP с сингулярной обработкой

ODE является единичным при r=0, поэтому нужно обрабатывать первый сегмент особым образом. Теорема о среднем значении дает 

(y'(r)-y'(0))/r->y''(0)  for  r->0,

, так что в этом пределе r->0 вы получите 

3*y''(0) = a*y(0) + b*y(0)^3`.

Это позволяет определить первый ar c как 

y(r) = y0 + (a*y0 + b*y0^3)*r^2/6
y'(r) = (a*y0 + b*y0^3)*r/3

до заказа r³ и r². Таким образом, если вам нужна точность 1e-9 в y(r), то первый сегмент не должен быть длиннее 1e-3.

Вместо того, чтобы пытаться исключить y0 из уравнений для y(h) и y'(h) чтобы получить условие, связывающее ya[0] и ya[1], пусть решатель также выполнит эту работу и добавит y0 в качестве параметра в систему. Тогда граничные условия имеют 3 слота, соответствующих добавленному виртуальному измерению для параметра, который может быть естественным образом заполнен уравнениями y(h)=ya[0] и ya[1]=y'(h) и правым граничным условием.

В целом вы можете определить систему как 

h = 1e-3;

def fun(r, y, p):
    return  y[1], -2/r*y[1]+y[0]*(a+b*y[0]*y[0]) 

def bc(ya, yb, p):
    y0, = p
    yh = y0 + y0*(a+b*y0*y0)*h*h/6;
    dyh = y0*(a+b*y0*y0)*h/3
    return ya[0]-yh, ya[1]-dyh, yb[0]-c


x = np.linspace(h, 10, 10)
ya = np.zeros((2, x.size))

sol = solve_bvp(fun, bc, x, ya, p=[1])
print(sol.message,f"y(0)={sol.p[0]}");
plt.plot(sol.x, sol.y[0]);

, так что с примерами параметров a, b, c = -1, 0.2, 3 вы получите конвергентный вызов решателя с y(0)=2.236081087849196 и результирующий график

enter image description here

...