我这样做是为了使用射击方法解决 BVP 函数。
我不明白为什么当我的步长 h = 0.1 时，限制 x_end 超出了我在函数solve_ivp下的 while 循环中设置的值。
当我在 while 循环中设置 x_end = 1 和 h = 0.1 时
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()))
我在 while 循环内得到 x 的以下迭代：
0.1 0.2 0.30000000000000004 0.4 0.5 0.6 0.7 0.7999999999999999 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 步长，但没有看到任何错误。
0.1 0.2 0.30000000000000004 0.4 0.5 0.6 0.7 0.7999999999999999 0.8999999999999999 0.9999999999999999 1.0999999999999999
x 增加到 1.09999999999999999 因为
要解决#1，您可以避免使用浮点数或实现一个函数来确定两个数字是否足够接近：
def almost_equal(a, b):
return abs(a**2-b**2)<1e-6
对于 #2，您可以将所有逻辑移至递增 x 的行之前，以避免 x 值出现这种差异。
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)
solution.append((x, y.copy()))
print(x)
x += h