Я надеюсь использовать scipy.integrate.solve_bvp для решения дифференциального уравнения 2-го порядка: я проверяю свой процесс с помощью предыдущего уравнения, поэтому я уверен, что перейду к более сложным уравнениям.
Начнем с системы дифференциальных уравнений:
f''(x) + f(x) - f(x)^3 = 0
с учетом граничных условий
f(x=0) = 0 f(x->infty) = gammaA
где gammaA - некоторая константа между 0 и 1. Я нахожу для этого численные решения и сравниваю с известной аналитической формой (по крайней мере, для gammaA = 1, функцией tanh). Для любого заданного gammaA мы можем один раз проинтегрировать это уравнение и использовать BC на бесконечности.
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp
gammaA = 0.9
xstart = 0.0
xend = 10
steps = 0.1
x = np.arange(xstart,xend,steps)
def dpsidx3(x,psi, gammaA):
eq = ( gammaA**2 *(1 - (1/2)*gammaA**2) - psi**2 *(1 - (1/2)*psi**2) )**0.5
return eq
psi0 = 0
x0 = xstart
x1 = xend
sol = solve_ivp(dpsidx3, [x0, x1], y0 = [psi0], args = (gammaA,), dense_output=True, rtol = 1e-9)
plotsol = sol.sol(x)
plt.plot(x, plotsol.T,marker = "", linestyle="--",label = r"Numerical solution - $solve\_ivp$")
plt.xlabel('x')
plt.ylabel('psi')
plt.legend()
plt.show()
Если gammaA не равно 1, то есть некоторые предупреждения во время выполнения, но форма точно такая, как ожидалось. Однако ОДУ в коде solve_ivp был преобразован в форму, которая является ОДУ 1-го порядка; для дальнейшей работы (с более сложными и переменными коэффициентами в ОДУ) это будет невозможно. Следовательно, я пытаюсь решить проблему граничного значения, используя solve_bvp. Я пытаюсь решить сейчас ту же ОДУ, но не получаю такого же результата, как от этого решения; в документации неясно, как эффективно использовать solve_bvp для меня! Вот моя попытка:
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_bvp
gammaA = 0.9
xstart = 0.0
xend = 10
steps = 0.1
def fun(x,u):
du1 = u[1] #u[1] = u2, u[0] = u1
du2 = u[0]**3 - u[0]
return np.vstack( (du1, du2) )
def bc(ua, ub):
return np.array( [ua[0], ub[0]-gammaA])
x = np.linspace(xstart, xend, 10)
print(x.size)
y_a = np.zeros((2, x.size))
y_a[0] = np.linspace(0, gammaA, 10)
y_a[0] = gammaA
res_a = solve_bvp(fun, bc, x, y_a, max_nodes=100000, tol=1e-9)
print(res_a)
x_plot = np.linspace(0, xend, 100)
y_plot_a = res_a.sol(x_plot)[0]
fig2,ax2= plt.subplots()
ax2.plot(x_plot, y_plot_a, label=r'BVP solve')
ax2.legend()
ax2.set_xlabel("x")
ax2.set_ylabel("psi")
Я попытался написать ОДУ 2-го порядка как систему ОДУ 1-го порядка и установить правильные граничные условия в конце системы (а не в бесконечности). Я ожидал подобной tanh-функции (где я мог бы сказать, что после окончания системы мое решение будет просто gammaA, как и ожидалось по асимптоте), но ясно, что я не получаю этого ни при каком значении gammaA. Любые советы с благодарностью; как я могу воспроизвести результат solve_ivp в solve_bvp?
РЕДАКТИРОВАТЬ: дополнительные мысли. Могу ли я добавить к своей задаче дополнительное ограничение, чтобы гарантировать, что решение имеет стационарную точку на краю/является монотонно возрастающим решением? Графики выглядят нормально для gammaA =1, но не показывают правильного поведения для любых других значений, как в solve_ivp.
EDIT2: сравнительные данные, показывающие плохое соответствие с gammaA, например. 0,8, но хорошее совпадение для gammaA = 1. 
< img src="https://i.stack.imgur.com/GLJ90.png" alt="gamma1" />
