specialForElena

import numpy as np
import matplotlib.pyplot as plt

# Define the problem parameters
def p(x):
    return 1.0

def q(x):
    return 1.0

def r(x):
    return x

def f(x):
    return np.exp(x)

# Set up the grid
a = 0.0
b = 1.0
N = 10  # initial number of grid points
x = np.linspace(a, b, N)

# Discretize the differential equation
h = x[1] - x[0]
A = np.zeros((N, N))
b = np.zeros(N)
for i in range(1, N - 1):
    A[i, i - 1] = -p(x[i]) / h**2 - q(x[i]) / (2 * h)
    A[i, i] = 2 * p(x[i]) / h**2 + r(x[i])
    A[i, i + 1] = -p(x[i]) / h**2 + q(x[i]) / (2 * h)
    b[i] = f(x[i])

# Set up the boundary conditions
A[0, 0] = 1.0
A[N - 1, N - 1] = 1.0
b[0] = 0.0
b[N - 1] = 0.0

# Solve the linear system
y = np.linalg.solve(A, b)

# Refine the grid using Richardson extrapolation
M = 2 * N  # number of points on the refined grid
x_ = np.linspace(a, b, M)
h_ = x_[1] - x_[0]
A_ = np.zeros((M, M))
b_ = np.zeros(M)
for i in range(1, M - 1):
    A_[i, i - 1] = -p(x_[i]) / h_**2 - q(x_[i]) / (2 * h_)
    A_[i, i] = 2 * p(x_[i]) / h_**2 + r(x_[i])
    A_[i, i + 1] = -p(x_[i]) / h_**2 + q(x_[i]) / (2 * h_)
    b_[i] = f(x_[i])
A_[0, 0] = 1.0
A_[M - 1, M - 1] = 1.0
b_[0] = 0.0
b_[M - 1] = 0.0
y_ = np.linalg.solve(A_, b_)

# Check the convergence rate
err = np.abs(y_ - y[::2])
h = h_[::2]
plt.loglog(h, err, '-o')
plt.xlabel('h')
plt.ylabel('error')
plt.show()

# Output the solution
plt.plot(x, y, '-o', label='coarse grid')
plt.plot(x_, y_, '-o', label='refined grid')
plt.legend()
plt.xlabel('x')
plt.ylabel('y')
plt.show()