Итерация выходит за пределы цикла while ⇐ Python

[Anonymous]

Я сделал это для решения функции BVP с использованием метода съемки.
Я не понимаю, почему, когда мой размер шага h = 0,1, предел x_end выходит за пределы того, что я установил в Цикл while под функциейsolve_ivp.
Когда я установил x_end = 1 и h = 0,1 в цикле while

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

x = 0
y = np.array(y0)
solution = [(x, y.copy())]

while x < x_end:
y = runge_kutta_4(f, x, y, h, p, q, r)
x += h
print(x)
solution.append((x, y.copy()))

Я получаю следующие итерации x внутри цикла while:
0.1
0.2
0.30000000000000004
0.4
0,5
0,6
0,7
0,79999999999999999
0.8999999999999999
0.9999999999999999
1.0999999999999999
Другие значения h не вызывают такой проблемы.
Может кто-нибудь объяснить почему и как я могу это исправить?
Ниже приведен мой полный код, если это поможет.

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

import numpy as np
import matplotlib.pyplot as plt

# Define the ODE as a system of first-order equations
def ode_system(x, y, p, q, r):
"""
Converts the second-order ODE into a system of two first-order ODEs.
y = [y1, y2], where y1 = y and y2 = y'.
"""
y1, y2 = y
dy1_dx = y2
dy2_dx = -p(x) * y2 - q(x) * y1 + r(x)
return np.array([dy1_dx, dy2_dx])

def runge_kutta_4(f, x0, y0, h, p, q, r):
"""
Fourth-order Runge-Kutta method for solving the system of ODEs.
"""
k1 = h * f(x0, y0, p, q, r)
k2 = h * f(x0 + h / 2, y0 + k1 / 2, p, q, r)
k3 = h * f(x0 + h / 2, y0 + k2 / 2, p, q, r)
k4 = h * f(x0 + h, y0 + k3, p, q, r)
return y0 + (k1 + 2 * k2 + 2 * k3 + k4) / 6

def solve_ivp(f, x_range, y0, h, p, q, r):
"""
Solves the initial value problem using the RK4 method.
"""
x0, x_end = x_range
x = x0
y = np.array(y0)
solution = [(x, y.copy())]

while x < x_end:
y = runge_kutta_4(f, x, y, h, p, q, r)
x += h
solution.append((x, y.copy()))

return solution

def shooting_method(a, b, alpha1, beta1, gamma1, alpha2, beta2, gamma2, p, q, r, h, tol=1e-6):
"""
Solves the boundary value problem using the shooting method.
- a, b: Interval [a, b]
- alpha1, beta1, gamma1: Boundary condition at x=a (third kind)
- alpha2, beta2, gamma2: Boundary condition at x=b (third kind)
- p, q, r: Coefficients in the ODE
- h: Step size
- tol: Tolerance for root finding
"""

def boundary_residual(guess):
y_a = (gamma1 - alpha1 * guess)/beta1
y0 = [y_a, guess]  # Initial conditions [y(a), y'(a)]
solution = solve_ivp(ode_system, (a, b), y0, h, p, q, r)
y_b, y_b_prime = solution[-1][1]  # Values of y(b) and y'(b)
return alpha2 * y_b_prime + beta2 * y_b - gamma2

# Initial guesses for y'(a)
g1, g2 = 0.5, 1.0  # Two initial guesses for the secant method

while abs(g2 - g1) > tol:
r1 = boundary_residual(g1)
r2 = boundary_residual(g2)
if abs(r2 - r1) <  tol:
raise ValueError("Secant method failed: residuals too close.")

# Secant update
g_new = g2 - r2 * (g2 - g1) / (r2 - r1)
g1, g2 = g2, g_new

# Solve the IVP with the best guess
y_best = [(gamma1 - alpha1 * g2)/beta1, g2]

return solve_ivp(ode_system, (a, b), y_best, h, p, q, r)

# Define the ODE: y'' + p(x)y' + q(x)y = r(x)
# Example: y"(x) - e**x * y'(x) - x * y(x) = 0
# Coefficients for the ODE
def p(x): return -np.exp(x)
def q(x): return -x
def r(x): return 0

# Third-kind boundary conditions
# alpha1 * y'(a) + beta1 * y(a) = gamma1
# alpha2 * y'(b) + beta2 * y(b) = gamma2
a, b = 0, 1 # x limits
alpha1, beta1, gamma1 = 1, -1, 1  # Boundary condition at x=0: y'(0) - y(0)  = 1
alpha2, beta2, gamma2 = 1, 1, 0  # Boundary condition at x=1: y'(1) + y(1) = 0

# Solve the BVP
h = 0.1  # Step size
solution = shooting_method(a, b, alpha1, beta1, gamma1, alpha2, beta2, gamma2, p, q, r, h)

# Extract the solution
x_vals = [x for x, y in solution]
y_vals = [y[0] for x, y in solution]

# Plot the solution
plt.plot(x_vals, y_vals, label="y(x)")
plt.xlabel("x")
plt.ylabel("y(x)")
plt.title("Solution of the BVP using Shooting Method")
plt.legend()
plt.grid()
plt.show()

Я попробовал изменить размер шага h, ошибок не обнаружено.

Подробнее здесь: https://stackoverflow.com/questions/792 ... while-loop