使用Python solve_bvp求解边值问题时遇除法无效值错误求助

Hey there, let's break down why you're hitting that division error and how to fix it—your problem is definitely tied to the M and dM equations you added, so let's focus there.

First, pinpoint the error source #

Looking at your fun function, the d2M calculation has multiple terms dividing by T:

d2M = -(dM/T)*dT + (dM/T)*theta*(m+1) - (alpha/T)*B*dB

When the solver iterates, if T drops to 0 or becomes NaN, you'll get that division error. And since removing the M equations fixes things, this is the clear culprit.

Top fixes to try #

1. Fix your initial guess (this is the biggest issue!) #

Right now, you're initializing all y values to 0:

y = np.zeros((8, x.size))

That means your initial T (which is y[4]) starts at 0—so you're dividing by 0 immediately! No wonder the solver throws an error.

Instead, set initial values that match your boundary conditions:

• T needs to go from 1 (left boundary) to 4 (right boundary)
• M needs to go from 0 (left) to 1 (right)
• For dM, start with a reasonable slope (since M goes from 0 to 1 over a 1-unit interval, a slope of 1 makes sense)

Update your initial guess code like this:

x = np.linspace(-0.5, 0.5, 500)
y = np.zeros((8, x.size))
# Initialize T with a linear transition from 1 to 4
y[4] = np.linspace(1, 4, x.size)
# Initialize M with a linear transition from 0 to 1
y[6] = np.linspace(0, 1, x.size)
# Initialize dM with a constant slope (matches the linear M guess)
y[7] = np.full(x.size, 1.0)

2. Add a safety guard against T approaching 0 #

Even with a good initial guess, the solver might push T close to 0 during iteration. Add a tiny epsilon to T to avoid division by zero:

def fun(x, y):
    U, dU, B, dB, T, dT, M, dM = y;
    # Small epsilon to prevent division by zero/near-zero T
    eps = 1e-8
    T_safe = T + eps
    d2U = -2*U_0*Q**2*(1/np.cosh(Q*x))**2*np.tanh(Q*x)-((alpha)/(C_k*sigma))*dB;
    d2B = -(1/(C_k*zeta))*dU;
    d2T = (1/gamma - 1)*(sigma*dU**2 + zeta*alpha*dB**2);
    # Use T_safe instead of T in all division terms
    d2M = -(dM/T_safe)*dT + (dM/T_safe)*theta*(m+1) - (alpha/T_safe)*B*dB
    # Make sure to return a numpy array (helps with solver stability)
    return np.array([dU, d2U, dB, d2B, dT, d2T, dM, d2M])

Stick with an epsilon like 1e-8—small enough that it doesn't skew your results, but large enough to avoid numerical issues.

3. Verify your boundary conditions #

Quick check: your boundary condition function returns a list, but it's better to return a numpy array for consistency. Also, double-check that each condition is correctly set to 0 (since solve_bvp expects residual values to be zero):

def bc(ya, yb):
    return np.array([
        ya[0] + U_0*np.tanh(Q*0.5),  # Left U boundary: ya[0] = -U0*tanh(Q*0.5)
        yb[0] - U_0*np.tanh(Q*0.5),  # Right U boundary: yb[0] = U0*tanh(Q*0.5)
        ya[2] - 0,                   # Left B boundary: ya[2] = 0
        yb[2] - 0,                   # Right B boundary: yb[2] = 0
        ya[4] - 1,                   # Left T boundary: ya[4] = 1
        yb[4] - 4,                   # Right T boundary: yb[4] = 4
        ya[6],                       # Left M boundary: ya[6] = 0
        yb[6] - 1                    # Right M boundary: yb[6] = 1
    ])

4. Tweak solver parameters if needed #

If you still run into issues, adjust the solver's settings to give it more flexibility:

sol = solve_bvp(fun, bc, x, y, max_nodes=1000, tol=1e-6)

• max_nodes lets the solver use more points to refine the solution
• tol adjusts the convergence tolerance (smaller = more precise, but slower)

Full working code example #

Here's the complete adjusted code with all fixes:

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

alpha = 1
zeta = 1
C_k = 1
sigma = 1
Q = 30
U_0 = 0.1
gamma = 5/3
theta = 3
m = 1.5

def fun(x, y):
    U, dU, B, dB, T, dT, M, dM = y;
    eps = 1e-8
    T_safe = T + eps
    d2U = -2*U_0*Q**2*(1/np.cosh(Q*x))**2*np.tanh(Q*x)-((alpha)/(C_k*sigma))*dB;
    d2B = -(1/(C_k*zeta))*dU;
    d2T = (1/gamma - 1)*(sigma*dU**2 + zeta*alpha*dB**2);
    d2M = -(dM/T_safe)*dT + (dM/T_safe)*theta*(m+1) - (alpha/T_safe)*B*dB
    return np.array([dU, d2U, dB, d2B, dT, d2T, dM, d2M])

def bc(ya, yb):
    return np.array([
        ya[0] + U_0*np.tanh(Q*0.5),
        yb[0] - U_0*np.tanh(Q*0.5),
        ya[2] - 0,
        yb[2] - 0,
        ya[4] - 1,
        yb[4] - 4,
        ya[6],
        yb[6] - 1
    ])

x = np.linspace(-0.5, 0.5, 500)
y = np.zeros((8, x.size))
# Set reasonable initial guesses
y[4] = np.linspace(1, 4, x.size)
y[6] = np.linspace(0, 1, x.size)
y[7] = np.full(x.size, 1.0)

sol = solve_bvp(fun, bc, x, y)

# Check if solver succeeded and plot results
if sol.success:
    plt.figure(figsize=(12, 8))
    plt.subplot(221)
    plt.plot(sol.x, sol.y[0], label='U')
    plt.title('U(x)')
    plt.legend()
    plt.subplot(222)
    plt.plot(sol.x, sol.y[2], label='B')
    plt.title('B(x)')
    plt.legend()
    plt.subplot(223)
    plt.plot(sol.x, sol.y[4], label='T')
    plt.title('T(x)')
    plt.legend()
    plt.subplot(224)
    plt.plot(sol.x, sol.y[6], label='M')
    plt.title('M(x)')
    plt.legend()
    plt.tight_layout()
    plt.show()
else:
    print(f"Solver failed: {sol.message}")

Quick debugging tips #

• Always check sol.success and sol.message to see why the solver failed
• Plot your initial guess first to make sure it aligns with boundary conditions
• For stiff systems or sensitive equations, good initial guesses are half the battle

内容的提问来源于stack exchange，提问作者Patrick Lewis