求解耦合二阶 ODE，Python 中的数值

Question

我想在 python 中对以下 DGL 系统进行数值求解：

程序应与往常相同。我有 2 个二阶耦合微分方程，我使用替换 g' = v 和 f' = u 创建四个一阶耦合微分方程。

这是我的代码：

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

# Constants
e = 1
mu = 1
lambd = 1

# Initial second-order system

# g'' + 2/r g'-2/r^2 g - 3/r eg^2- e^2 g^3 - e/(2r) f^2 - e^2/4 f^2g = 0
# f'' + 2/r f' - 2/r^2 f - 2e/r fg - e^2/2 fg^2 + mu^2 f - lambda f^3 = 0

# Substitution
# g' = v
# v' = g''
# f' = u
# u' = f''

# Equivalent set of first-order systems
def dSdr(r, S):
    g, v, f, u = S
    dgdr = v
    dvdr = -2 / r * v + 2 / r ** 2 * g + 3 / r * e * g ** 2 + \
        e ** 2 * g ** 3 + e / (2 * r) * f**2 + e**2 / 4 * f ** 2 * g
    dfdr = u
    dudr = -2 / r * u + 2 / r ** 2 * f + 2 * e / r * f * g + e ** 2 / 2 * f * g**2 - mu ** 2 * f + lambd * f ** 3

    return [dgdr, dvdr, dfdr, dudr]

# Initial conditions
g0 = 0.0001
v0 = 0.0001
f0 = 1
u0 = 0.0001
S0 = [g0, v0, f0, u0]

r_start = 0.1
r_end = 100
r = np.linspace(r_start, r_end, 1000)

# Solve the differential equations

sol = solve_ivp(dSdr, t_span=(min(r), max(r)), y0=S0, t_eval=r, method='RK45')




# Check if integration was successful
if sol.success:
    print("Integration successful")
else:
    print("Integration failed:", sol.message)

# Debugging information
print("Number of function evaluations:", sol.nfev)
print("Number of accepted steps:", sol.t_events)
print("Number of steps that led to errors:", sol.t)

# Plot results
plt.figure(figsize=(10, 6))
plt.plot(sol.t, sol.y[0], label='g(r)')
plt.plot(sol.t, sol.y[2], label='f(r)')
plt.xlabel('r')
plt.ylabel('Function values')
plt.legend()
plt.title('Solutions of the differential equations')
plt.show()

我的边界骗局。应该是 f(0) = g(0) = 0 且 f(infinity) = mu/sqrt(lambd) g(infinity) = 0 以便系统具有物理意义。但是我如何合并这个条件或者我如何知道 v,u 的初始条件？系统应该看起来像这样（来自原始论文）：

但事实并非如此。有谁知道该怎么办吗？

来源：

在此输入链接描述

Answer 1

虽然您可以尝试使用初始值求解器的“射击”方法，但我发现最好直接解决边界值问题。我通过类似于流体力学中的有限体积方法来实现这一点。

除以 r 会产生奇点，因此首先将方程乘以 r²。然后您将得到看起来像球对称扩散方程的方程：

这里，S_f和S_g是很棒的长源术语（您可能需要在我的代码中检查它们 - 有很多小错误的范围）。将每个离散化（相当于通过有限体积方法对单元 [r_i-1/2,r_i+1/2] 进行积分）。例如对于 f：

这个以及 g 的相应方程可以写成以下形式

边界处存在微小差异，f 上的非零边界条件产生了非平凡的解。

然后可以通过三对角矩阵算法或高斯-赛德尔或雅可比等迭代方案迭代求解离散方程。为了简单起见，我在这里选择了稍微松弛的雅可比方法。

import numpy as np
import matplotlib.pyplot as plt

N = 100                                # number of cells
rmin, rmax = 0.0, 20.0                 # minimum and maximum r
dr = rmax / N                          # cell width
gL, gR, fL, fR = 0, 0, 0, 1            # boundary conditions

e, nu, lamda = 1, 1, 1

r = np.linspace( rmin + 0.5 * dr, rmax - 0.5 * dr, N )       # cell centres
aw = np.zeros( N )
ap = np.zeros( N )
ae = np.zeros( N )
sg = np.zeros( N )
sf = np.zeros( N )
sp = np.zeros( N )

# Connection coefficients at boundaries
awL = 2 * rmin ** 2 / dr ** 2
aeR = 2 * rmax ** 2 / dr ** 2

for i in range( N ):
    sp[i] = -2
    if i == 0:
        aw[i] = 0
        sp[i] -= awL
    else:
        aw[i] = ( 0.5 * ( r[i-1] + r[i] ) ) ** 2 / dr ** 2

    if i == N - 1:
        ae[i] = 0
        sp[i] -= aeR
    else:
        ae[i] = ( 0.5 * ( r[i+1] + r[i] ) ) ** 2 / dr ** 2
    ap[i] = aw[i] + ae[i] - sp[i]

# Initial values
g = np.zeros( N )
f = np.zeros( N )
gtarget = g.copy()
ftarget = f.copy()
f = f + 1
g = g + 1

alpha = 0.9                  # Under-relaxation factor

niter = 0
for _ in range( 10000 ):
    niter += 1

    # Standard source term
    sg = -3*e*r*g**2 -     e**2*r**2*g**3  -0.5*e*r*f**2 - 0.25*e**2*r**2*f**2*g
    sf = -2*e*r*f*g  - 0.5*e**2*r**2*f*g**2+    nu**2*r**2*f-lamda*r**2*f**3

    # Dirichlet boundary conditions applied via source term
    sg[0  ] += awL * gL
    sf[0  ] += awL * fL
    sg[N-1] += aeR * gR
    sf[N-1] += aeR * fR

    # Update g and f (under-relaxed Jacobi method)
    gtarget = g.copy()
    ftarget = f.copy()
    for i in range( 1, N - 1 ):
        gtarget[i] = ( sg[i] + aw[i] * g[i-1] + ae[i] * g[i+1] ) / ap[i]
        ftarget[i] = ( sf[i] + aw[i] * f[i-1] + ae[i] * f[i+1] ) / ap[i]
    i = 0
    gtarget[i] = ( sg[i] + ae[i] * g[i+1] ) / ap[i]
    ftarget[i] = ( sf[i] + ae[i] * f[i+1] ) / ap[i]
    i = N - 1
    gtarget[i] = ( sg[i] + aw[i] * g[i-1] ) / ap[i]
    ftarget[i] = ( sf[i] + aw[i] * f[i-1] ) / ap[i]
    g = alpha * gtarget + ( 1 - alpha ) * g
    f = alpha * ftarget + ( 1 - alpha ) * f
    error = np.linalg.norm( f - ftarget ) + np.linalg.norm( g - gtarget )
    # print( niter, error )
    if error < 1.0e-6: break

print( niter, " iterations " )

plt.plot( r, f, label='f' )
plt.plot( r, g, label='g' )
plt.grid()
plt.legend()
plt.show()