如何使用 python 数值求解具有多个边界条件的二阶导数方程

Question

正如标题所建议的，现在我正在尝试使用 python 数值求解二阶导数方程。使用的Python包是Scipy，特别是solve_bvp。我正在使用的方程是

d²y/dx²-ay-by³=0，其中 < 0 and b > 0.

有四个边界条件：y(±∞) = ±1 和 dy/dx = 0，在 x=±∞ 处，根据 本文，方程的精确解为 y ∝ tanh(x) 形式可以获得。所以我想做的是得到与这个精确解类似的数值结果。

下面提供了我通过chatgpt编写和修改的代码。使用solve_bvp，我能够包含四个边界条件中的两个，即y(±∞) = ±1。然而，我得到的数值结果根本不像 tanh 函数，而是像 ，我怀疑这是因为我在使用solve_bvp 时只能提供 2 个边界条件。当我尝试包含四个边界条件时，系统总是发出一条错误消息：

ValueError: `bc` return is expected to have shape (2,), but actually has (4,).

我想知道是否有一种方法可以将额外的两个边界条件添加到solve_bvp中，或者有其他求解器可以支持多个边界条件。

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

a = -1
b = 1
p0 = 1

# Define the differential equation
def y_derivative(x, y):
    return np.vstack((y[1], a * y[0] + b * (y[0]**3)))

# Define the boundary conditions
def bc(ya, yb):
    return np.array([ya[0] + p0, yb[0] - p0])

# Generate a set of points
x = np.linspace(-100, 100, 100)

# Solve the boundary value problem
y0 = np.zeros((2, x.size))
sol = solve_bvp(y_derivative, bc, x, y0)

# Calculate the exact solution
y_exact = np.tanh(x)

# Plot the numerical solution
plt.plot(x, sol.sol(x)[0], label='Numerical Solution')

# Plot the exact solution
plt.plot(x, y_exact, label='Exact Solution', linestyle='--')

# Customize plot
plt.xlabel('x')
plt.ylabel('y')
plt.title('Comparison of Numerical and Exact Solutions')
plt.legend()
plt.grid(True)
plt.show()

Answer 1

让我们以您的 MCVE 为基础：

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

注意我们添加了参数

p

:

def model(x, y, p):
    return np.array([
        y[1],
        p[0] * y[0] + p[1] * y[0] ** 3
    ])

我们还需要将其添加到

bc

：

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

我们对域名进行了一些限制：

x = np.linspace(-5., 5., 500)

并准备网格使其更接近真实情况：

y0 = np.zeros((2, x.size))
y0[0, :200] = -1.
y0[0, -200:] = +1.

现在集成表现良好：

sol = integrate.solve_bvp(model, bc, x, y0, p=(-1., 1.), max_nodes=5000)

#       message: The algorithm converged to the desired accuracy.
#       success: True
#        status: 0
#             x: [-5.000e+00 -4.980e+00 ...  4.980e+00  5.000e+00]
#           sol: <scipy.interpolate._interpolate.PPoly object at 0x7f66e55ff740>
#             p: [-8.642e+00  8.642e+00]
#             y: [[-1.000e+00 -1.000e+00 ...  1.000e+00  1.000e+00]
#                 [ 0.000e+00  1.299e-09 ...  1.299e-09  0.000e+00]]
#            yp: [[ 0.000e+00  1.299e-09 ...  1.299e-09  0.000e+00]
#                 [ 6.476e-08  6.499e-08 ... -6.499e-08 -6.476e-08]]
# rms_residuals: [ 5.282e-14  5.863e-14 ...  5.978e-14  5.293e-14]
#         niter: 2