用solve_bvp检查solve_ivp的结果 - solve_bvp问题

编程问答 2022-05-02

问题描述

我希望使用 scipy.integrate.solve_bvp 来求解二阶微分方程：我正在用之前的方程检查我的过程，所以我有信心继续研究更复杂的方程。

我们从微分方程组开始：

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 代码中的 ODE 已被处理为一阶 ODE 形式；对于进一步的工作（在 ODE 中使用更复杂和可变的系数），这是不可能的。因此，我正在尝试使用 solve_bvp 解决边界值问题。 我现在正在尝试解决相同的 ODE，但我没有得到与此解决方案相同的结果；文档不清楚如何有效地使用 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,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,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")

我尝试将二阶 ODE 编写为一阶 ODE 系统，并在系统末端（而不是无穷远）设置正确的边界条件。我期望有一个类似的 tanh 函数（我可以说在系统结束后，我的解决方案只是 gammaA，正如渐近线所预期的那样），但很明显我没有得到任何价值gammaA。非常感谢任何建议；如何在 solve_ivp 中重现 solve_bvp 的结果？

编辑：额外的想法。 我可以为我的问题添加额外的约束以确保解决方案在边缘有一个固定点/是一个单调递增的解决方案吗？对于 gammaA =1，这些图看起来不错，但对于 solve_ivp 中的任何其他值都没有显示正确的行为。

EDIT2：比较数字，显示与 gammaA 的一致性较差，例如0.8 但对于 gammaA = 1 的一致性很好。

gamma0.8

gamma1

解决方法

您对这个方程的数学性质做出了毫无根据的假设。有一个能量泛函

E = u'^2 + u^2 - 0.5*u^4 - 0.5 = u'^2 - 0.5*(u^2-1)^2

您首先计算的解是在 0 能级。

对于任何较小的负能级，大致在单位圆内，您会得到周期性振荡解，这些解在无穷远处没有限制。 对于更大的正能级，解是无界的，将迅速变大，可能在有限的时间内发散。同样在这里，无穷大的极限要么不存在，因为没有将初始点连接到大时间的解决方案，要么极限就是无穷大本身。

强制边界条件违反这种性质可能会起作用，但不会给出稳定的解决方案。

boundarydifferential-equationsodepythonscipyscipy