Scipy 两点边界值问题

计算科学 scipy
2021-12-02 04:59:53

非线性 ODE 语句 我想使用 scipy 来解决以下问题:

u'' + (u')^2 = sin(x)
u(0)=0,
u(1)=1    

其中 u = u(x)。

方法 我正在查看以下链接

https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.solve_bvp.html

问题 - 不确定变量 wrt x 或如何输入 BC

链接中给出的示例不包含 wrt x 一词,特别是不确定如何包含

sin(x)

从我的问题的右侧,也不是

u(1)=1

完成工作

在下面,我尝试修改上面链接中的示例以符合我上面的 ODE

   from scipy.integrate import *
import numpy as np

from numpy import *
def fun(x, y):
    return np.vstack((y[1],np.sin(x) -(y[0])**2))
def bc(ya, yb):
    return np.array([ya[0], yb[0]-1])
x = np.linspace(0, 1, 5)
y_a = np.zeros((2, x.size))
y_b = np.zeros((2, x.size))
y_b[0] = 3

from scipy.integrate import solve_bvp
res_a = solve_bvp(fun, bc, x, y_a)
res_b = solve_bvp(fun, bc, x, y_b)
x_plot = np.linspace(0, 1, 100)
y_plot_a = res_a.sol(x_plot)[0]
y_plot_b = res_b.sol(x_plot)[0]
import matplotlib.pyplot as plt
plt.plot(x_plot, y_plot_a, label='y_a')
plt.plot(x_plot, y_plot_b, label='y_b')
plt.legend()
plt.xlabel("x")
plt.ylabel("y")
plt.show()
1个回答

我认为在 NumPy 中,如果y是一个矩阵,y[0]它会为您提供第 0 行,但如果您使用y[0,:].

我看到的主要问题是fun你在定义 ODE

u+u2=sinx,
使用y[0,:]而不是y[1,:].