Решение нелинейной системы ОДУ с граничными условиями на Python ⇐ Python

[Anonymous]

Я пытаюсь решить эту систему дифференциальных уравнений:


Функции должны вести себя так, как r-> бесконечность (F = g = лямбда = const.):
< img alt="введите описание изображения здесь" src="https://i.sstatic.net/AWVDL68J.png" />
Некоторое время назад у меня возникла аналогичная проблема, и пользователь был настолько любезен, что предоставил мне алгоритм, с помощью которого я мог бы решить эту ODS. Теперь я попытался применить алгоритм TDM к новой проблеме, но получил только тривиальное решение. Кто-нибудь знает, что я сделал не так или как я могу решить эту проблему?
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
p = -c / denominator
q = d / denominator
for i in range(1, n):
denominator = b + a * p
if abs(denominator) < 1.0e-10:
print("No solution")
return x
p = -c / 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()


Подробнее здесь: https://stackoverflow.com/questions/790 ... -in-python