具有奇异性和无穷大边界条件的边值问题

Question

我试图解决0，\infty上的以下边值问题

F‘=-f’/r+f/r^2+m^2*f+2 \lambda *f^3

f(0)=0 \；f(\infty)=\sqrt{-m^2/(2\lambda)}

对于一些常量m^2<0，\lambda>0。不存在闭型，但f应该从0单调增加到sqrt{-m^2/(2\lambda)。在r=0中有一个可移动的奇点。这个问题就是贝塞尔方程加上f^3中的一个项。

我试图用integrate.solve_bvp来解决这个问题，它可以解决在一个边界上具有奇异性的多边界问题，定义y=f，rf‘，以便

y'=0，r(m^2f+2\lambda f^3)+(1/r)*[0,1,1,0]*y

将边界条件设为无穷大值max_x。不幸的是，我的代码遵循了https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.solve_bvp.html的示例结构，给出了错误的解决方案：

import scipy.integrate
import numpy as np
import matplotlib.pyplot as plt

m_squared=-1
Lambda=1
asymptote=np.sqrt(-m_squared/(2*Lambda))

#evaluate infinity b.c here
max_x=100

def fun(r,v):
    z=(m_squared*v[0]+(2*Lambda)*(v[0]**3) )*r
    return np.vstack((z-z, z))
  
#boundary condition
def bc(ya,yb):
    return np.array([ya[0], yb[0]-asymptote])

# to treat singularity
S=np.array([[0,1],[1,0]])


x=np.linspace(0,max_x,5000)

# guess for vector y at points x  
y=np.zeros((2, len(x)))
y[0,-1]=asymptote
print(y)

#solve
res=scipy.integrate.solve_bvp(fun, bc, x, y, p=None, S=S)
x_plot=np.linspace(0,max_x,1000)
y_plot=res.sol(x_plot)[0]

plt.plot(x_plot,y_plot,label="numerical")
plt.axhline(asymptote,linestyle="--",label="asymptote")

plt.xlabel("r")
plt.ylabel("f")
plt.legend()

​

​

我检查了修改上面的代码以解决$f''=f-1$和$f(0)=0，f(\infty)=1$很好。在这种情况下不存在奇异性，因此修改fun和设置S=None就足够了。

我的代码是否有问题，还是应该使用不同的边值求解器？

Answer 1

我从你的来源中读到的方程式是

f''(r)+f'(r)/r-f(r)/r^2 = 2*lambda*eta^2*f(r)*(f(r)^2-1)

使用建议的参数lambda=0.2, eta=2。

在长期极限中，左侧降为二阶导数，方程化为一个守恒系统，中心在f=0，两个鞍点在f=+-1。任务是找到一条收敛到鞍点的解曲线。在更实际的条件下，这类似于推动一个刚性钟摆的方式，使它以直立的姿势结束，或者更接近于那个位置。

为接近f=1鞍点的解编写f=1，该方程近似为

g''(r) = a^2*g(r),  a^2=4*lambda*eta^2=3.2

这再次将这个平衡描述为一个鞍点，向它收敛的解满足约化的ODE g'(r)=-a*g(r)。这可以作为上边界条件。转换成状态向量，这给出

def bc(ya,yb):
    return np.array([ya[0], yb[1]+a*max_x*(yb[0]-1)])

(如果您想保留您的版本，请用asymptote替换平衡常数1)。

我从这方面得到了很好的结果，包括论文中的参数和你的参数。

然而，求解器的奇异模式似乎被打破了，它将节点插入到允许的max_nodes接近于零的位置，其中的解应该是简单线性的。我把最初的猜测设为

x=np.logspace(-5,np.log10(max_x),10)
x[0]=0

# guess for vector y at points x  
y=[np.tanh(a*x),a*x/np.cosh(a*x)**2]

这样就不会从一开始就违反这个最大节点数。

不使用奇异机制，将奇异项转化为ODE函数，利用有限解在初始段几乎是线性的情况下，可以使用f(r)=r*f'(r)作为初始条件，ya[0]-ya[1] == 0。然后区间开始是一些小的正数。这将导致解决方案中合理的节点数目，默认公差为1e-3为25，公差为1e-6为100。