Python求解带边界条件的非线性ODE方程组遇平凡解问题

我尝试求解一组非线性常微分方程组，当r趋近于无穷时，函数需满足如下渐近行为：

F(r) → 1，W(r) → 1/r

此前我用TDMA算法成功解决过类似问题，但本次应用该算法仅得到平凡解，恳请协助排查错误或提供可行的解决方法，以下是我的Python实现代码：

import numpy as np
import matplotlib.pyplot as plt

# TDMA function to solve the tri-diagonal matrix equation
def tdma(a, b, c, d):
    n = len(d)
    p = np.zeros(n)
    q = np.zeros(n)
    x = np.zeros(n)

    # Forward pass
    i = 0
    denominator = b[i]
    p[i] = -c[i] / denominator
    q[i] = d[i] / denominator
    for i in range(1, n):
        denominator = b[i] + a[i] * p[i - 1]
        if abs(denominator) < 1.0e-10:
            print("No solution")
            return x
        p[i] = -c[i] / denominator
        q[i] = (d[i] - a[i] * q[i - 1]) / denominator

    # Backward pass
    x[n - 1] = q[n - 1]
    for i in range(n - 2, -1, -1):
        x[i] = p[i] * x[i + 1] + q[i]

    return x


# Parameters
N = 1000  # Number of cells (1-N); added boundaries 0 and N+1
rmin, rmax = 0.0, 100.0  # Minimum and maximum r
dr = rmax / N  # Cell width
lamda = 10
g = 1.0  # You can adjust this parameter if needed

# Boundary conditions on F and W
FL, FR = 0, 1
WL, WR = 0, 1 / rmax

# Cell / node arrangement
r = np.linspace(rmin - 0.5 * dr, rmax + 0.5 * dr, N + 2)
r[0], r[N + 1] = rmin, rmax
rv = np.linspace(rmin, rmax, N + 1)

# Coefficients
awF = np.zeros(N + 2)
aeF = np.zeros(N + 2)
awW = np.zeros(N + 2)
aeW = np.zeros(N + 2)

# Main cells
aeF[1:N + 1] = 1 / dr ** 2
aeW[1:N + 1] = rv[1:N + 1] ** 2 / dr ** 2
awF[1:N + 1] = aeF[0:N]
awW[1:N + 1] = aeW[0:N]

# Boundaries
awF[1], awW[1], aeF[N], aeW[N] = 0, 0, 0, 0
awFL = 2 / dr ** 2
aeFR = 2 / dr ** 2
awWL = 2 * rmin ** 2 / dr ** 2
aeWR = 2 * rmax ** 2 / dr ** 2

# Initial values
F = np.ones(N + 2) * 1e-5
W = np.ones(N + 2) * 1e-5

# Boundary conditions
F[0], F[N + 1] = FL, FR
W[0], W[N + 1] = WL, WR

# Under-relaxation factor
alpha = 0.5
niter = 0
s = np.zeros(N + 2)
sp = np.zeros(N + 2)

# Iteration loop
for _ in range(2000):
    niter += 1
    F_target = F.copy()
    W_target = W.copy()

    ######################
    # Update W equation  #
    ######################
    for i in range(1, N + 1):
        K = 1 - g * r[i] * W[i]
        H = g * r[i] * F[i]

        # Clip values to prevent overflow
        K = np.clip(K, -1e5, 1e5)
        H = np.clip(H, -1e5, 1e5)

        s[i] = K * (K**2 - 1) + H**2 * K
        sp[i] = -2 * (1 - g * r[i]) / r[i]**2

        if np.isnan(K) or np.isnan(H):
            print(f"Numerical issue at r = {r[i]}, K = {K}, H = {H}")
            break

    s[1] += awFL * FL
    sp[1] -= awFL
    s[N] += aeFR * FR
    sp[N] -= aeFR

    ap = awF + aeF - sp
    F_target[1:N + 1] = tdma(-awF[1:N + 1], ap[1:N + 1], -aeF[1:N + 1], s[1:N + 1])

    ######################
    # Update F equation  #
    ######################
    for i in range(1, N + 1):
        K = 1 - g * r[i] * W[i]
        H = g * r[i] * F[i]

        # Clip values to prevent overflow
        K = np.clip(K, -1e5, 1e5)
        H = np.clip(H, -1e5, 1e5)

        s[i] = 2 * H * K**2 + lamda * H * ((H**2 / g**2) - r[i]**2 * F[i]**2)

        if np.isnan(K) or np.isnan(H):
            print(f"Numerical issue at r = {r[i]}, K = {K}, H = {H}")
            break

    s[1] += awWL * WL
    sp[1] -= awWL
    s[N] += aeWR * WR
    sp[N] -= aeWR

    ap = awW + aeW - sp
    W_target[1:N + 1] = tdma(-awW[1:N + 1], ap[1:N + 1], -aeW[1:N + 1], s[1:N + 1])

    # Calculate change
    change = (np.linalg.norm(F_target - F) + np.linalg.norm(W_target - W)) / N

    # Under-relax the update
    F = alpha * F_target + (1 - alpha) * F
    W = alpha * W_target + (1 - alpha) * W

    if change < 1.0e-10:
        break

print(niter, " iterations ")

# Plot the results
plt.plot(r, F, label='F(r)', color='g')
plt.plot(r, W, label='W(r)', color='r')
plt.legend(loc="upper right")
plt.xlabel("r")
plt.ylim(0)
plt.xlim(0, 20)
plt.grid(True)
plt.show()

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