Approximate u''(x) = sin^2(x/2), 0 < x < 2*Pi with boundary conditions u(0) = 0 and u(2*Pi) = 2*Pi. using: (1) 2nd-order accurate finite-difference approximation; (2) 4th-order accurate finite-difference approximation; (3) 4th-order extrapolation method on two grids; (4) 4th-order deferred correction; - spectral collocation Plot the error as h, …

fd_bvp.py

import fd.fd_stencil as fd
import numpy as np
import math
import matplotlib.pyplot as plt

"""
Approximate u''(x) = sin^2(x/2), 0 < x < 2*Pi
with boundary conditions u(0) = 0 and u(2*Pi) = 2*Pi. 
using: 
        - 2nd-order accurate finite-difference approximation;
        - 4th-order accurate finite-difference approximation;
        - 4th-order extrapolation method on two grids;
        - 4th-order deferred correction;
        - spectral collocation
Plot the error as h, the mesh width, goes to zero. 
"""

# setting up the mesh width
number_entries = 18
h_values = np.logspace(0, -3, number_entries)
interval_upper = 2 * math.pi


# setting up the function
def f(x):
    return math.sin(x/2)**2


def u(x):
    return 1/4 * (x**2 + 2 * math.cos(x))


# (1) solving the 2nd-order FD approximation
fd_coefficients = fd.fd_stencil(2, [-1, 0, 1])
errors = np.zeros(number_entries)

counter = 0
for h in h_values:
    n_points = math.ceil(interval_upper/h)
    n_intervals = n_points - 1
    actual_h = interval_upper / n_intervals

    n_vars = n_points - 2  # two boundary conditions

    x_values = np.array([i * actual_h for i in range(0, n_points)])
    solution = np.array([u(x_i) for x_i in x_values])

    matrix_A = np.zeros((n_vars, n_vars))
    rhs = np.zeros(n_vars)

    matrix_A[0, 0:2] = fd_coefficients[1:3]  # first row
    matrix_A[n_vars-1, n_vars-2:n_vars] = fd_coefficients[0:2]  # last row
    rhs[0] = f(1.0 * actual_h) - 0.5/actual_h**2
    rhs[n_vars-1] = f(n_vars * actual_h) - (0.5 + math.pi**2)/actual_h**2

    for i in range(1, n_vars - 1):  # excludes first and last rows
        matrix_A[i, i-1:i+2] = fd_coefficients
        x_i = (i + 1) * actual_h
        rhs[i] = f(x_i)

    matrix_A = 1/actual_h**2 * matrix_A
    answer = np.linalg.solve(matrix_A, rhs)
    answer = np.insert(answer, 0, 0.5)
    answer = np.insert(answer, n_vars + 1, 0.5 + math.pi**2)

    calculated_error = np.abs(answer - solution)
    errors[counter] = np.linalg.norm(calculated_error, np.inf)
    counter = counter + 1


# (2) solving the 4th-order FD approximation
fd_coefficients = fd.fd_stencil(2, [-2, -1, 0, 1, 2])
errors2 = np.zeros(number_entries)

counter = 0
for h in h_values:
    n_points = math.ceil(interval_upper/h)
    n_intervals = n_points - 1
    actual_h = interval_upper / n_intervals

    n_vars = n_points - 2  # two boundary conditions

    x_values = np.array([i * actual_h for i in range(0, n_points)])
    solution = np.array([u(x_i) for x_i in x_values])

    matrix_A = np.zeros((n_vars, n_vars))
    rhs = np.zeros(n_vars)

    matrix_A[0, 0:5] = fd.fd_stencil(2, [-1, 0, 1, 2, 3, 4])[1:]  # first row
    matrix_A[1, 0:4] = fd.fd_stencil(2, [-2, -1, 0, 1, 2])[1:]  # second row
    matrix_A[n_vars-2, n_vars-4:n_vars] = fd.fd_stencil(2, [-2, -1, 0, 1, 2])[:4]  # second-to-last row
    matrix_A[n_vars-1, n_vars-5:n_vars] = fd.fd_coefficients(2, 0, [-4, -3, -2, -1, 0, 1])[:5]  # last row

    rhs[0] = f(1.0 * actual_h) - (10 * 0.5)/(12.0 * actual_h**2)
    rhs[1] = f(1.0 * actual_h) + 0.5 / (12.0 * actual_h**2)
    rhs[n_vars-2] = f((n_vars-1) * actual_h) + (0.5 + math.pi**2) / (12.0 * actual_h**2)
    rhs[n_vars-1] = f(n_vars * actual_h) - (10 * (0.5 + math.pi ** 2)) / (12.0 * actual_h**2)

    for i in range(2, n_vars - 2):  # excludes two firsts and two last rows
        matrix_A[i, i-2:i+3] = fd_coefficients
        x_i = (i + 1) * actual_h
        rhs[i] = f(x_i)

    matrix_A = (1.0 / actual_h**2) * matrix_A
    answer = np.linalg.solve(matrix_A, rhs)
    answer = np.insert(answer, 0, 0.5)
    answer = np.insert(answer, n_vars + 1, 0.5 + math.pi**2)

    calculated_error = np.abs(answer - solution)
    errors2[counter] = np.linalg.norm(calculated_error, np.inf)
    counter = counter + 1


# (3) solving the 4th-order extrapolation method
fd_coefficients = fd.fd_stencil(2, [-1, 0, 1])
errors3 = np.zeros(number_entries)
# to check validity of each individual solution
errors_fine = np.zeros(number_entries)
errors_coarse = np.zeros(number_entries)
h_values_coarse = 2.0 * h_values

counter = 0
for h in h_values:
    # coarse grid
    h_coarse = h_values_coarse[counter]
    n_points_coarse = math.ceil(interval_upper / h_coarse)
    n_intervals_coarse = n_points_coarse - 1
    actual_h_coarse = interval_upper / n_intervals_coarse
    n_vars_coarse = n_points_coarse - 2  # two boundary conditions

    matrix_A_coarse = np.zeros((n_vars_coarse, n_vars_coarse))
    rhs_coarse = np.zeros(n_vars_coarse)

    matrix_A_coarse[0, 0:2] = fd_coefficients[1:3]  # first row
    matrix_A_coarse[n_vars_coarse - 1, n_vars_coarse - 2:n_vars_coarse] = fd_coefficients[0:2]  # last row
    rhs_coarse[0] = f(1.0 * actual_h_coarse) - 0.5 / actual_h_coarse ** 2
    rhs_coarse[n_vars_coarse - 1] = f(n_vars_coarse * actual_h_coarse) - (0.5 + math.pi ** 2) / actual_h_coarse ** 2

    for i in range(1, n_vars_coarse - 1):  # excludes first and last rows
        matrix_A_coarse[i, i - 1:i + 2] = fd_coefficients
        x_i = (i + 1) * actual_h_coarse
        rhs_coarse[i] = f(x_i)

    matrix_A_coarse = (1 / actual_h_coarse ** 2) * matrix_A_coarse
    answer_coarse = np.linalg.solve(matrix_A_coarse, rhs_coarse)

    # check coarse grid solution for accuracy
    x_values_coarse = np.array([i * actual_h_coarse for i in range(0, n_points_coarse)])
    solution_coarse = np.array([u(x_i) for x_i in x_values_coarse])
    answer_coarse = np.insert(answer_coarse, 0, 0.5)
    answer_coarse = np.insert(answer_coarse, n_vars_coarse + 1, 0.5 + math.pi ** 2)
    calculated_error_coarse = np.abs(answer_coarse - solution_coarse)
    errors_coarse[counter] = np.linalg.norm(calculated_error_coarse, np.inf)

    # fine grid
    n_points = 2 * n_points_coarse - 1
    n_intervals = n_points - 1
    actual_h = interval_upper / n_intervals

    n_vars = n_points - 2  # two boundary conditions

    matrix_A = np.zeros((n_vars, n_vars))
    rhs = np.zeros(n_vars)

    matrix_A[0, 0:2] = fd_coefficients[1:3]  # first row
    matrix_A[n_vars - 1, n_vars - 2:n_vars] = fd_coefficients[0:2]  # last row
    rhs[0] = f(1.0 * actual_h) - 0.5 / actual_h ** 2
    rhs[n_vars - 1] = f(n_vars * actual_h) - (0.5 + math.pi ** 2) / actual_h ** 2

    for i in range(1, n_vars - 1):  # excludes first and last rows
        matrix_A[i, i - 1:i + 2] = fd_coefficients
        x_i = (i + 1) * actual_h
        rhs[i] = f(x_i)

    matrix_A = (1 / actual_h**2) * matrix_A
    answer_fine = np.linalg.solve(matrix_A, rhs)

    x_values = np.array([i * actual_h for i in range(0, n_points)])
    solution = np.array([u(x_i) for x_i in x_values])

    # test for fine grid accuracy
    answer_fine = np.insert(answer_fine, 0, 0.5)
    answer_fine = np.insert(answer_fine, n_vars + 1, 0.5 + math.pi**2)
    calculated_error = np.abs(answer_fine - solution)
    errors_fine[counter] = np.linalg.norm(calculated_error, np.inf)

    # combine the fine and coarse solutions
    combined_answer = np.zeros(answer_coarse.size)

    for i in range(0, answer_coarse.size):
        combined_answer[i] = 1.0/3.0 * (4.0 * answer_fine[2 * i] - answer_coarse[i])

    # for i in range(0, answer_fine.size):
    #     combined_answer[i] = 1/3 * (4 * answer_fine[i] - answer_coarse[min(math.ceil(i/2), answer_coarse.size-1)])

    # combined_answer = np.insert(combined_answer, 0, 0.5)
    # combined_answer = np.insert(combined_answer, combined_answer.size, 0.5 + math.pi ** 2)

    calculated_error = np.abs(combined_answer - solution_coarse)
    errors3[counter] = np.linalg.norm(calculated_error, np.inf)
    counter = counter + 1
    print(errors3)


# (4) solving the 4th-order deferred correction method
fd_coefficients = fd.fd_stencil(2, [-1, 0, 1])
errors4 = np.zeros(number_entries)


def d2f(x):
    return math.cos(x) / 2.0


def correction(x, h):
    return 1/12.0 * h**2 * d2f(x)


counter = 0
for h in h_values:
    n_points = math.ceil(interval_upper/h)
    n_intervals = n_points - 1
    actual_h = interval_upper / n_intervals

    n_vars = n_points - 2  # two boundary conditions

    x_values = np.array([i * actual_h for i in range(0, n_points)])
    solution = np.array([u(x_i) for x_i in x_values])

    matrix_A = np.zeros((n_vars, n_vars))
    rhs = np.zeros(n_vars)

    matrix_A[0, 0:2] = fd_coefficients[1:3]  # first row
    matrix_A[n_vars-1, n_vars-2:n_vars] = fd_coefficients[0:2]  # last row
    rhs[0] = f(1.0 * actual_h) - 0.5/actual_h**2 + correction(1.0 * actual_h, actual_h)
    rhs[n_vars-1] = f(n_vars * actual_h) - (0.5 + math.pi**2)/actual_h**2 + correction(n_vars * actual_h, actual_h)

    for i in range(1, n_vars - 1):  # excludes first and last rows
        matrix_A[i, i-1:i+2] = fd_coefficients
        x_i = (i + 1) * actual_h
        rhs[i] = f(x_i) + correction(x_i, actual_h)

    matrix_A = 1/actual_h**2 * matrix_A
    answer = np.linalg.solve(matrix_A, rhs)
    answer = np.insert(answer, 0, 0.5)
    answer = np.insert(answer, n_vars + 1, 0.5 + math.pi**2)

    calculated_error = np.abs(answer - solution)
    errors4[counter] = np.linalg.norm(calculated_error, np.inf)
    counter = counter + 1


# (5) solving pseudospectral method
def get_chebyshev_point(a, b, z):  # x = ax + (bx-ax) * 0.5*(1 + cos(pi*(1-z)));
    return a + (b - a) * 0.5 * (1 + math.cos(math.pi * (1 - z)))


counter = 0
num_points = [8, 12, 16, 22]
errors5 = np.zeros(len(num_points))
ps_h_values = np.zeros(len(num_points))
for n_points in num_points:
    actual_h = interval_upper / (n_points - 1)

    matrix_A = np.zeros((n_points+1, n_points+1))
    rhs = np.zeros(n_points+1)

    x_entries = np.linspace(0, 1, n_points + 1)  # x_entries = np.array([i for i in range(0, n_points-1)])
    x_values = np.array([get_chebyshev_point(0, 2 * math.pi, x) for x in x_entries])
    matrix_A[0, :] = fd.fd_coefficients(0, x_values[0], x_values)  # first row
    matrix_A[n_points, :] = fd.fd_coefficients(0, x_values[n_points], x_values)  # last row

    rhs[0] = 0.5
    rhs[n_points] = 0.5 + math.pi**2

    for i in range(1, n_points):
        matrix_A[i, :] = fd.fd_coefficients(2, x_values[i], x_values)
        rhs[i] = f(x_values[i])

    answer = np.linalg.solve(matrix_A, rhs)

    solution = np.array([u(x_i) for x_i in x_values])
    calculated_error = np.abs(answer - solution)
    errors5[counter] = np.linalg.norm(calculated_error, np.inf)
    ps_h_values[counter] = actual_h
    counter = counter + 1

print(errors5)

# plot
plt.figure(figsize=(8, 6))
plt.title("Order of convergence for FD methods")
plt.loglog(h_values, errors, alpha=0.35, marker='o', label="2nd-order FD")
plt.loglog(h_values, errors2, alpha=0.65, marker='s', label="4th-order FD")
plt.loglog(h_values, errors3, alpha=0.35, marker='x', label="Extrapolation")
plt.loglog(h_values, errors4, alpha=0.5, marker='2', label="Deferred correction")
plt.loglog(ps_h_values, errors5, alpha=0.5, marker='D', label="Pseudospectral")
# to check individual solutions for the extrapolation method
# plt.loglog(h_values, errors_fine, alpha=0.35, marker='.', label="Extrapolation-Fine")
# plt.loglog(h_values_coarse, errors_coarse, alpha=0.35, marker='o', label="Extrapolation-Coarse")

plt.xlabel("Mesh width (h)")
plt.ylabel("Approximation error")
plt.grid()
plt.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()