python - 用solve_bvp检查solve_ivp的结果——solve_bvp问题
问题描述
我希望用它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, 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")
我试图将二阶 ODE 编写为一阶 ODE 的系统,并在系统末尾(而不是无穷大)设置正确的边界条件。我期望有一个类似的 tanh 函数(我可以说在系统结束后,我的解决方案很简单gammaA
,正如渐近线所预期的那样),但很明显,对于gammaA
. 任何建议都将不胜感激;如何重现solve_ivp
in的结果solve_bvp
?
编辑:额外的想法。我可以为我的问题添加一个额外的约束以确保解决方案在边缘有一个固定点/是一个单调递增的解决方案吗?这些图看起来不错gammaA =1
,但没有显示任何其他值的正确行为,如solve_ivp
.
解决方案
你对这个方程的数学性质做出了毫无根据的假设。存在能量泛函
E = u'^2 + u^2 - 0.5*u^4 - 0.5 = u'^2 - 0.5*(u^2-1)^2
您首先计算的解决方案位于能级 0。
对于任何较小的负能级,大致在单位圆内,你会得到周期性的振荡解,它们在无穷大没有限制。对于更大的正能量水平,解决方案是无界的,将迅速变大,可能在有限时间内发散。同样在这里,无穷大的极限要么不存在,因为没有解决方案将初始点与大时间联系起来,要么极限本身就是无穷大。
针对这种性质强制边界条件可能会起作用,但不会给出稳定的解决方案。
推荐阅读
- javascript - 用 Node Js 初始化 PostgreSQL
- ios - Alamofire.upload - 未调用响应处理程序
- php - Codecept Laravel 模拟 mysql 数据库将值转换为 json
- javascript - 如何从 jquery 对话框中添加和删除表
- angular - 需要同时禁用复选框并使其看起来被禁用(样式))
- laravel - 在 LARAVEL 5.5 中的单个操作中将表列更新为 1,将其他行更新为 0
- haproxy - 如何在 haproxy 中创建粘性会话
- razor - 谁能解释一下 .split(''); 在剃刀?
- kubernetes - 在 pod 之间共享文件
- python - 在 Python 中使用另一个类的函数和数据