Как завершить метод стрельбы по проблеме с ОДУ? ⇐ Python

[Anonymous]

Сейчас я пытаюсь найти метод съемки с использованием Python.
Исходный вопрос:
ODE, с m=0
Вот мой код, и мне интересно почему график станет более плоским, если допуск уменьшится.

Код: Выделить всё

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

# index
B = 4.0
I = 20.65
sigma = 2.372
tolerance = 1e-10

# eqn
def ode_system(x, y, a):
f, fp = y
if x < 1e-8:  # avoid 0
fpp = -a * sigma
else:
r = x / (1 - x)
r2 = r**2
denom = 1 + 2 * B * np.sin(f)**2 / r2
num = -fp * (2 / r) - np.sin(2 * f) / r2 * (B * (fp**2 - 1) - I * np.sin(f)**2 / r2)
fpp = num / denom
fpp *= (1 - x)**4
fpp -= 2 * (1 - x)**3 * fp
return [fp, fpp]

# return f(1) , f'(1)
def boundary_shooting(a):
f0 = np.pi
fp0 = 0  #  f'(0) = 0
y0 = [f0, fp0]

sol = solve_ivp(
ode_system,
[0, 0.99999999],  # avoid x=1
y0,
args=(a,),
t_eval=[0.99999999],
rtol=tolerance,
atol=tolerance,
)
if not sol.success:
print("Integration failed:", sol.message)
f_at_1 = sol.y[0, -1]  # f(1)
fp_at_1 = sol.y[1, -1]  # f'(1)
return f_at_1, fp_at_1

# find a
def find_optimal_a():
a_low, a_high = 0.000000011769865256, 0.000000011769865257
iteration = 0
while a_high - a_low > tolerance:
a_mid = (a_low + a_high) / 2
f_end, fp_end = boundary_shooting(a_mid)

# show the message
print(f"Iteration {iteration}: a_low={a_low:.15e}, a_high={a_high:.15e}, a_mid={a_mid:.15e}, f(1)={f_end:.15e}, f'(1)={fp_end:.15e}")
iteration += 1

if abs(f_end) < tolerance and abs(fp_end) < tolerance:  # let f(1),f'(1)=0
return a_mid

if f_end > 0:
a_high = a_mid
else:
a_low = a_mid

return (a_low + a_high) / 2

# main
if __name__ == "__main__":
# find a
a_optimal = find_optimal_a()
print("Optimal a:", a_optimal)

f0 = np.pi
fp0 = 0  # make f'(0) = 0
y0 = [f0, fp0]
sol = solve_ivp(
ode_system,
[0, 0.99999999],
y0,
args=(a_optimal,),
t_eval=np.linspace(0, 0.99999999, 100000),
rtol=tolerance,
atol=tolerance,
)

# boudary value
x_vals = sol.t
f_vals = sol.y[0]
fp_vals = sol.y[1]

print("f(0):", f_vals[0])
print("f'(0):", fp_vals[0])
print("f(1):", f_vals[-1])
print("f'(1):", fp_vals[-1])

# Plot
plt.plot(x_vals, f_vals, label="f(x)")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.title("Shooting Method Solution in x-coordinates")
plt.legend()
plt.grid()
plt.show()

Я пытался изменить допуск, но не смог получить стабильный график функции f(x). Я хотел бы увидеть результат, что f(1)~0 и f'(1)~0. Я не совсем понимаю, где я ошибся, но самый нормальный результат, который я получил, был
f(0): 3.141592653589793
f'(0): 0.0f''(0): -2.7918120388418003e-08
f(1): 0,9460917810683763
f'(1): -2,231818513581139

Подробнее здесь: https://stackoverflow.com/questions/793 ... de-problem