import numpy as np
import matplotlib.pyplot as plt
def ode(x, y, dy):
“””Second order ODE to solve: y” = f(x,y,y’)”””
return 4*y – 4*x # Example ODE: y” = 4y – 4x
def rk4_step(f, x, y, dy, h):
“””Single step of RK4 for second order ODE”””
# First set of weights
k1 = h * dy
l1 = h * f(x, y, dy)
# Second set of weights
k2 = h * (dy + l1/2)
l2 = h * f(x + h/2, y + k1/2, dy + l1/2)
# Third set of weights
k3 = h * (dy + l2/2)
l3 = h * f(x + h/2, y + k2/2, dy + l2/2)
# Fourth set of weights
k4 = h * (dy + l3)
l4 = h * f(x + h, y + k3, dy + l3)
# Update y and dy
y_new = y + (k1 + 2*k2 + 2*k3 + k4)/6
dy_new = dy + (l1 + 2*l2 + 2*l3 + l4)/6
return y_new, dy_new
def shooting_method(f, x0, y0, x_end, y_end, dy_guess1, dy_guess2, h, tol=1e-6, max_iter=100):
“””Shooting method implementation”””
# First shot
x = np.arange(x0, x_end + h, h)
y1 = np.zeros_like(x)
dy1 = np.zeros_like(x)
y1[0] = y0
dy1[0] = dy_guess1
for i in range(len(x)-1):
y1[i+1], dy1[i+1] = rk4_step(f, x[i], y1[i], dy1[i], h)
# Second shot
y2 = np.zeros_like(x)
dy2 = np.zeros_like(x)
y2[0] = y0
dy2[0] = dy_guess2
for i in range(len(x)-1):
y2[i+1], dy2[i+1] = rk4_step(f, x[i], y2[i], dy2[i], h)
# Secant method iterations
for _ in range(max_iter):
# New guess using secant method
dy_new = dy_guess2 – (y2[-1] – y_end) * (dy_guess2 – dy_guess1) / (y2[-1] – y1[-1])
# Update guesses
dy_guess1, dy_guess2 = dy_guess2, dy_new
# Compute new trajectory
y1, y2 = y2, np.zeros_like(x)
dy2 = np.zeros_like(x)
y2[0] = y0
dy2[0] = dy_guess2
for i in range(len(x)-1):
y2[i+1], dy2[i+1] = rk4_step(f, x[i], y2[i], dy2[i], h)
# Check convergence
if abs(y2[-1] – y_end) < tol:
break
return x, y2
# Boundary conditions
x0, y0 = 0, 0 # y(0) = 0
x_end, y_end = 1, 1 # y(1) = 1
# Initial guesses for y'(0)
dy_guess1 = 0
dy_guess2 = 1
# Step size
h = 0.1
# Solve using shooting method
x, y = shooting_method(ode, x0, y0, x_end, y_end, dy_guess1, dy_guess2, h)
# Exact solution for comparison (for the example ODE)
exact = lambda x: x + np.sinh(2*x)/np.sinh(2)
y_exact = exact(x)
# Print results
print(f"{'x': < 8}{'Numerical': < 12}{'Exact': < 12}{'Error': < 12}")
print("-" * 40)
for xi, yi, ye in zip(x, y, y_exact):
print(f"{xi: < 8.3f}{yi: < 12.6f}{ye: < 12.6f}{abs(yi-ye): < 12.6f}")
# Plot results
plt.figure(figsize=(10, 6))
plt.plot(x, y, 'bo-', label='Numerical Solution')
plt.plot(x, y_exact, 'r--', label='Exact Solution')
plt.title("Solution of Boundary Value Problem using Shooting Method")
plt.xlabel("x", fontsize=12)
plt.ylabel("y(x)", fontsize=12)
plt.legend(fontsize=10, edgecolor="black")
plt.grid(True, linestyle="--")
plt.show()
