Naomi D.
Yamamoto
1/31/25
BSCPE - 3 Numerical
Methods
Problem 1: Fixed-Point Iteration Method
Problem:
The Fixed-Point Iteration Method is a numerical approach for finding roots of the
equation
g(x)=xg(x) = x. Your task is to implement the Fixed-Point Iteration Method and visualize
the
iterations.
Requirements:
1. Implement the Fixed-Point Iteration Method.
2. Plot g(x)g(x), the line y=xy = x, and the iterations.
3. Inputs:
- g(x)g(x): The iterative function.
- x0x_0: The initial guess.
- max_iter: Maximum iterations (default: 10).
- tol: Tolerance for convergence (default: 10−610^{-6}).
4. Output:
- A plot showing g(x)g(x), y=xy = x, and the iterations.
- Print the approximate fixed point and the number of iterations.
Code:
# Import required libraries
import numpy as np
import matplotlib.pyplot as plt
# Copy all the functions
def fixed_point_iteration(g, x0, max_iter=10, tol=1e-6):
x_values = [x0]
x_old = x0
for i in range(max_iter):
x_new = g(x_old)
x_values.append(x_new)
if abs(x_new - x_old) < tol:
return x_values, i + 1
x_old = x_new
return x_values, max_iter
def plot_iterations(g, x_values, x_range):
plt.figure(figsize=(10, 8))
x = np.linspace(x_range[0], x_range[1], 1000)
plt.plot(x, g(x), 'b-', label='g(x)')
plt.plot(x, x, 'r--', label='y=x')
for i in range(len(x_values)-1):
plt.plot([x_values[i], x_values[i]],
[x_values[i], x_values[i+1]], 'g-', linewidth=1)
plt.plot([x_values[i], x_values[i+1]],
[x_values[i+1], x_values[i+1]], 'g-',
linewidth=1)
plt.plot(x_values, x_values, 'ko', label='Iterations')
plt.grid(True)
plt.legend()
plt.title('Fixed-Point Iteration Method')
plt.xlabel('x')
plt.ylabel('y')
plt.axis('equal')
plt.show()
# Define the function
def example_function(x):
return np.sqrt(10-x)
# Run in separate cells:
# Cell 1: Run with example_function
x0 = 1.0
max_iter = 10
tol = 1e-6
x_range = [0, 4]
x_values, iterations = fixed_point_iteration(example_function,
x0, max_iter, tol)
plot_iterations(example_function, x_values, x_range)
print(f"Approximate fixed point: {x_values[-1]:.6f}")
print(f"Number of iterations: {iterations}")
print(f"Iteration values: {[f'{x:.6f}' for x in x_values]}")
# Cell 2: Try with cos(x)
def new_function(x):
return np.cos(x)
x0 = 0.5
x_range = [-1, 1]
x_values, iterations = fixed_point_iteration(new_function, x0,
max_iter, tol)
plot_iterations(new_function, x_values, x_range)
print(f"Approximate fixed point: {x_values[-1]:.6f}")
print(f"Number of iterations: {iterations}")
print(f"Iteration values: {[f'{x:.6f}' for x in x_values]}")
Problem 2: Euler's Method for Solving ODEs
Problem:
Implement Euler's Method to approximate the solution of the differential equation
y′=f(x,y)y' =
f(x, y) over a given interval. Plot the numerical solution.
Requirements:
1. Implement Euler's Method.
2. Inputs:
- f(x,y)f(x, y): The function representing the derivative.
- x0,y0x_0, y_0: Initial condition.
- hh: Step size.
- nn: Number of steps.
3. Output:
- A plot of the numerical solution y(x)y(x).
- Print the final approximate value of yy.
Code:
import numpy as np
import matplotlib.pyplot as plt
# Differential Equation
def f(x, y):
return x - y
# Parameters
x0, y0 = 0, 1 # Initial condition
h = 0.1 # Step size
n = 20 # Number of steps
def euler_method(f, x0, y0, h, n):
x = np.zeros(n+1)
y = np.zeros(n+1)
x[0] = x0
y[0] = y0
for i in range(1, n+1):
x[i] = x[i-1] + h
y[i] = y[i-1] + h * f(x[i-1], y[i-1])
return x, y
# Run Euler Method
x, y = euler_method(f, x0, y0, h, n)
# Plotting
plt.figure(figsize=(10, 6))
plt.plot(x, y, 'b-', label="Numerical Solution")
plt.title("Euler's Method: Solution of y' = x - y")
plt.xlabel('x')
plt.ylabel('y')
plt.legend()
plt.grid(True)
plt.show()
# Print the final approximate value of y
print(f"Final approximate value of y: {y[-1]:.6f}")
Final approximate value of y: 1.243153
Problem 3: Simpson's Rule for Numerical Integration
Problem:
Implement Simpson's Rule to approximate the integral of a function f(x)f(x) over a given
interval. Visualize the function and the parabolas used in the approximation.
Requirements:
1. Implement Simpson's Rule.
2. Inputs:
- f(x)f(x): The function to integrate.
- a,ba, b: The integration limits.
- nn: Number of subintervals (must be even).
3. Output:
- A plot of f(x)f(x) with the parabolas.
- Print the approximate integral.
Code:
import matplotlib.pyplot as plt
import numpy as np
def simpsons_rule(f, a, b, n):
# Ensure n is even for Simpson's Rule
if n % 2 != 0:
raise ValueError("n must be even for Simpson's rule.")
# Variables
delta_x = (b - a) / n
x_values = np.linspace(a, b, n + 1) # Generate n+1 equally
spaced points
y_values = [f(x) for x in x_values] # Compute f(x) for each
point
# Compute Simpson's Rule
integral = f(a) + f(b) # First and last terms
for i in range(1, n):
coefficient = 4 if i % 2 != 0 else 2 # Alternate
coefficients
integral += coefficient * f(x_values[i])
integral *= (delta_x / 3)
# Print the result
print("Integral result:", integral)
# Plotting
x_plot = np.linspace(a, b, 500) # Smooth curve for f(x)
y_plot = f(x_plot)
# Plot the function
plt.plot(x_plot, y_plot, label="f(x) = x^3", color="blue")
# Highlight the integration points and area
plt.bar(x_values, y_values, width=delta_x, align='edge',
alpha=0.4, edgecolor="black", label="Integration Strips")
# Add labels and legend
plt.title("Simpson's Rule Visualization")
plt.xlabel("x")
plt.ylabel("f(x)")
plt.legend()
plt.grid(True)
# Show the plot
plt.show()
# Define the function f(x) = x^3
def cubic_function(x):
return x ** 3
# Run the program with f(x) = x^3
simpsons_rule(cubic_function, 0, 2, 10)