Решение двух связанных краевых задач второго порядка

Question

Я решил одно дифференциальное уравнение второго порядка с двумя граничными условиями, используя модуль solve_bvp. Однако сейчас я пытаюсь решить систему двух дифференциальных уравнений второго порядка:

U '' + a * B '= 0

B' '+ b * U' = 0

с граничными условиями U (+/- 0,5) = +/- 0,01 и B (+/- 0,5) = 0. Я разбил это на систему первых обыкновенных дифференциальных уравнений и я пытаюсь использовать solve_bvp для их численного решения. Тем не менее, я просто получаю массивы, полные нулей для моего решения. Я считаю, что неправильно выполняю граничные условия. Мне не понятно, как обращаться с более чем двумя уравнениями из документации. Моя попытка ниже

import numpy as np
from scipy.integrate import solve_bvp
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import solve_bvp
alpha = 1E-8
zeta = 8E-3
C_k = 0.05
sigma = 0.01
def fun(x, y):

    return np.vstack((y[1],-((alpha)/(C_k*sigma))*y[2],y[2], -(1/(C_k*zeta))*y[1]))

def bc(ya, yb):

    return np.array([ya[0]+0.001, yb[0]-0.001,ya[0]-0, yb[0]-0])

x = np.linspace(-0.5, 0.5, 5000)
y = np.zeros((4, x.size))
print(y)


sol = solve_bvp(fun, bc, x, y)
print(sol)

В моем вопросе я только что поменял a и b, но это просто параметры, которые я ввел. У меня есть аналитическое c решение для этой системы уравнений, поэтому я знаю, что существует не тривиальное. Любая помощь будет принята с благодарностью.

Answer 1

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

В форме вектор возврата производной, я думаю, вы намереваетесь

U, U',  B, B'

, что означает, что U=y[0], U'=y[1] и B=y[2],B'=y[3], так что ваш вектор производных должен быть правильно

return y[1], -((alpha)/(C_k*sigma))*y[3],  y[3], -(1/(C_k*zeta))*y[1]

и граница условия

return ya[0]+0.001, yb[0]-0.001, ya[2]-0, yb[2]-0

Особенно ваше граничное условие должно бросить алгоритм на первом шаге из-за особого якобиана, всегда проверять поле .success и поле .message структуры решения.

---

Обратите внимание, что по умолчанию абсолютный и относительный допуск экспериментального solve_bvp равен 1e-3, а количество узлов ограничено 500.

Установка начального номера узла равным 50 (5000 - это слишком много, решатель уточняет, где это необходимо) и допуск к 1-6, я получаю следующие графики решений, которые визуально удовлетворяют граничным условиям.