Решение связанного ОДУ 2-го порядка, числовое в Python
Я хотел бы численно решить следующую систему DGL в Python:
Процедура должна быть такой же, как всегда. У меня есть два связанных дифференциальных уравнения 2-го порядка, и я использую замену g' = v и f' = u, чтобы создать четыре связанных дифференциальных уравнения 1-го порядка.
Вот мой код:
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? Система должна выглядеть так (из оригинальной статьи):
но это не так. Кто-нибудь знает что делать?
Источник:
Я приложил источник, а также фрагмент статьи, где вы можете увидеть, какую систему они пытаются решить. Они масштабируют систему на коэффициент, прежде чем решить ее, но это не должно иметь значения, не так ли?
Ответы 1
Хотя вы можете попробовать метод «стрельбы» с использованием решателей начального значения, я считаю, что лучше сразу решить краевую задачу. Я делаю это методом, похожим на метод конечного объема в вычислительной механике жидкости (CFD).
Деление на r создает сингулярности, поэтому начните с умножения уравнений на r2. Собрав производные члены, вы получите уравнения, которые выглядят как сферически-симметричные уравнения диффузии:
Здесь Sf и Sg — отличные термины с длинным исходным кодом (которые вам, вероятно, придется проверить в моем коде — там много места для мелких ошибок). Дискретизируйте каждый (эквивалентно интегрированию по ячейке [ri-1/2,ri+1/2] методом конечного объема). например для ж:
Это и соответствующее уравнение для g можно записать в виде
На границах есть небольшие различия, и именно ненулевое граничное условие на f дает нетривиальное решение.
Дискретизированные уравнения затем могут быть решены итеративно либо с помощью алгоритма трехдиагональной матрицы (TDMA), либо с помощью итерационных схем, таких как Гаусс-Зейдель или Якоби. Хотя опыт работы с CFD подсказывает, что TDMA, вероятно, будет лучшим вариантом, для простоты я выбрал здесь (немного недостаточно расслабленный) метод Якоби.
Обратите внимание, что, как и в практике CFD, я в своем отредактированном коде разделил исходные термины, например: sf + spf.f, причем spf имеет отрицательное значение, поэтому его можно объединить с диагональным коэффициентом ap на LHS.
РЕДАКТИРОВАНИЕ: код изменен так, что значения центра ячейки равны 1...N с граничными индексами 0 и N+1, исходные термины разделены на явные и неявные части, код больше похож на Python, чем на Фортран!
import numpy as np
import matplotlib.pyplot as plt
N = 100 # number of cells (1-N); added boundaries 0 and N+1
rmin, rmax = 0.0, 20.0 # minimum and maximum r (latter an approximation to "infinity"
dr = rmax / N # cell width
gL, gR, fL, fR = 0, 0, 0, 1 # boundary conditions
e, nu, lamda = 1, 1, 1
# Cell / node arrangement:
#
# 0 N+1
# | 1 2 N-1 N |
# |-----|-----| . . . . . . |-----|-----|
#
r = np.linspace( rmin - 0.5 * dr, rmax + 0.5 * dr, N + 2 ) # cell centres (except for boundaries)
r[0] = rmin; r[N+1] = rmax # boundary nodes on faces of cells
# Connection coefficients
awL = 2 * rmin ** 2 / dr ** 2
aeR = 2 * rmax ** 2 / dr ** 2
aw = np.zeros( N + 2 )
ae = np.zeros( N + 2 )
for i in range( 1, N + 1 ):
if i == 1:
aw[i] = 0
else:
aw[i] = ( 0.5 * ( r[i-1] + r[i] ) ) ** 2 / dr ** 2
if i == N:
ae[i] = 0
else:
ae[i] = ( 0.5 * ( r[i+1] + r[i] ) ) ** 2 / dr ** 2
ap = aw + ae
# Initial values
g = np.zeros( N + 2 )
f = np.zeros( N + 2 )
f = f + 1
g = g + 1
# Boundary conditions
f[0] = fL; f[N+1] = fR; g[0] = gL; g[N+1] = gR
alpha = 0.9 # Under-relaxation factor
niter = 0
for _ in range( 10000 ):
niter += 1
# Source term (e.g., for f this would be decomposed as sf + spf.f, where spf is negative)
spf = -2 -0.5 * e**2 * r**2 * g**2 - lamda * r**2 * f**2
sf = -2 * e * r * f * g + nu**2 * r**2 * f
spg = -2 - e**2 * r**2 * g**2 - 0.25 * e**2 * r**2 * f**2
sg = -e * r * ( 3 * g**2 + 0.5 * f**2 )
# Dirichlet boundary conditions applied via source term
sf[1] += awL * fL; spf[1] -= awL
sg[1] += awL * gL; spg[1] -= awL
sf[N] += aeR * fR; spf[N] -= aeR
sg[N] += aeR * gR; spg[N] -= aeR
# Update g and f (under-relaxed Jacobi method)
ftarget = f.copy()
gtarget = g.copy()
for i in range( 2, N ): # cells 2 - (N-1)
ftarget[i] = ( sf[i] + aw[i] * f[i-1] + ae[i] * f[i+1] ) / ( ap[i] - spf[i] )
gtarget[i] = ( sg[i] + aw[i] * g[i-1] + ae[i] * g[i+1] ) / ( ap[i] - spg[i] )
i = 1 # leftmost cell
ftarget[i] = ( sf[i] + ae[i] * f[i+1] ) / ( ap[i] - spf[i] )
gtarget[i] = ( sg[i] + ae[i] * g[i+1] ) / ( ap[i] - spg[i] )
i = N # rightmost cell
ftarget[i] = ( sf[i] + aw[i] * f[i-1] ) / ( ap[i] - spf[i] )
gtarget[i] = ( sg[i] + aw[i] * g[i-1] ) / ( ap[i] - spg[i] )
# Under-relax the update
f = alpha * ftarget + ( 1 - alpha ) * f
g = alpha * gtarget + ( 1 - alpha ) * g
change = ( np.linalg.norm( f - ftarget ) + np.linalg.norm( g - gtarget ) ) / N
# print( niter, change )
if change < 1.0e-10: break
print( niter, " iterations " )
plt.plot( r, f, label='f' )
plt.plot( r, g, label='g' )
plt.grid()
plt.legend()
plt.show()
@Hendriksdf5, ваши уравнения в том виде, в котором вы их написали, везде полностью соответствуют тривиальному решению f=g=0. Так что, возможно, они ошибаются. Заявленные вами граничные условия также не соответствуют изображению, которое вы опубликовали. Пожалуйста, разместите ссылку на оригинальные уравнения, а не на вашу их версию. Откуда бы вы ни сделали эту фотографию, это будет началом. Уравнения такого типа можно решать методами конечного объема, но сначала нам нужно увидеть правильные уравнения.