我正在努力解决以下二阶边值问题:
y'' + 2/x*y' + k**2.0*F(y) = 0
y(x=1)=1, y'(x=0)=0
F(y) = -y or F(y) = -y*exp(AB*(1-y)/(1+B(1-y))
我以某种方式未能正确设置边界条件。我通过以下方式为F(y)=y
和边界条件定义了函数:
def fun(x, y, p):
k = p[0]
return np.vstack((y[1], -2.0/x*y[1] + k**2.0*y[0]))
def bc(ya, yb, p):
return np.array([ya[0], yb[0],ya[1]])
y[0,:] = 1
y[1,0] = 0
sol = solve_bvp(fun, bc, x, y, p=[40])
我应该得到的结果肯定是错误的,改变初始条件并不能使事情变得更好。我认为我的问题是与x = 0处的零梯度边界条件有什么关系。有人知道我在做什么错吗?
编辑:此处的MWE,对于k = 0.01,应给定常数1。但对于k = 5,x = 0处的值应约为1。 0.06:
def fun(x, y, p):
k = p[0]
return np.vstack((y[1], -2.0/x*y[1] + k**2.0*y[0]))
def bc(ya, yb, p):
return np.array([ya[0], yb[0]-1.0,yb[1]])
x = np.linspace(1e-3, 1, 100)
y = np.zeros((2, x.size))
y[0,:] = 1
from scipy.integrate import solve_bvp
sol = solve_bvp(fun, bc, x, y, p=[1000])
x_plot = np.linspace(0, 1, 100)
y_plot = sol.sol(x_plot)[0]
plt.figure()
plt.plot(x_plot, y_plot)
考虑情况F(y)=y
。这样就很容易看出该线性ODE的基本解是sin(k*x)/x
和cos(kx/x)
。这意味着大多数解在x=0
处具有奇异性。对于标准ODE求解器来说,开箱即用几乎是不可能的。必须在奇异点上帮助求解器,在离奇异点一定距离处再次应用常规过程。
您可以做的是分析x=0
处的情况,然后稍微走一下。您通过差异商的限制得到
y''(0) + 2*y''(0) + k^2*F(y(0)) = 0
允许您计算二次泰勒多项式。因此,使用到[a, 1]
上奇点y(x)=y(0)-k**2/6*F(y(0))*x**2
的延续,用ODE求解器解决[0,a]
上的问题。
x=a
的边界条件最容易建立,将y0=y(0)
作为参数。那么ODE和BC函数是
def ode(x,y,y0): return [ y[1], -2*y[1]/x - k**2*F(y[0]) ]
def bc(ya,yb,y0): y2 = -k**2*F(y0)/6; return [ ya[0] - y0 - y2*a**2, ya[1] - 2*y2*a, yb[0]-1 ]
在问题中讨论的情况下,这给出
a = 1e-2
def F(y): return y
for k in [0.01, 5]:
res = solve_bvp(ode, bc, [a,1], [[1,1], [0,0]], p=[1], tol=1e-5)
print(f'k={k}: {res.message}, y0={res.p[0]}, theory: {k/np.sin(k)}')
if res.success:
y0 = res.p[0]
x = np.linspace(a,1,61);
plt.plot(x, res.sol(x)[0])
plt.plot([0], [y0],'o', res.x, res.y[0],'+', ms=4)
plt.title(f'k={k}'); plt.grid(); plt.show()
结果
k = 0.01:算法收敛到所需的精度。,y0 = 1.00001666686111,理论值:1.0000166668611132
k = 5:该算法收敛到所需的精度。,y0 = -5.214175394024083,理论值:-5.214176063857029