gpt4 book ai didi

python - 使用 scipysolve_bvp 求解边值问题(扩散 react 方程)

转载 作者:太空宇宙 更新时间:2023-11-03 19:49:58 44 4
gpt4 key购买 nike

我正在努力解决以下二阶边界值问题:

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 处的值应约为。 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的情况。那么很容易看出这个线性常微分方程的基解是sin(k*x)/xcos(kx/x)。类似地,对于F(y)=-y,可以得到sinh(k*x)/xcosh(k*x)/x。这意味着大多数解在 x=0 处都有奇点。对于标准 ODE 求解器来说,这种奇点几乎不可能立即处理。人们必须在奇点处帮助求解器,在距奇点一定距离处再次应用正常程序。

你能做的就是分析x=0处的情况,然后稍微移开一点。你可以通过差商极限得到

y''(0) + 2*y''(0) + k^2*F(y(0)) = 0

它允许您计算二次泰勒多项式。因此,使用 ODE 求解器使用奇点延拓来求解 [a, 1] 上的问题 y(x)=y(0)-k**2/6*F(y (0))**x**2[0,a] 上。

y0=y(0) 视为参数,最容易建立x=a 处的边界条件。 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.sinh(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: The algorithm converged to the desired accuracy., y0=0.9999833335277726, theory: 0.9999833335277757

plot of solution

k=5: The algorithm converged to the desired accuracy., y0=0.06738256929427147, theory: 0.06738252915294544

plot of solution

关于python - 使用 scipysolve_bvp 求解边值问题(扩散 react 方程),我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/59904633/

44 4 0
Copyright 2021 - 2024 cfsdn All Rights Reserved 蜀ICP备2022000587号
广告合作:1813099741@qq.com 6ren.com