#226 Finite difference method to solve ODE using python

Solve the following problem:
d^2y/dt^2=-4*y + 4*x

Boundary conditions are
y(0) = 0 and y'(π/2) = 0.

The exact solution of the problem is
y = x − sin(2x*)

y(0) = 0,
y(i−1) − 2*y(i) + y(i+1) − h^2*(−4*y(i) + 4*x(i)) = 0, i = 1,2,…,n − 1
2*y(n−1) − 2*y(n) − h^2*(−4*y(n) + 4*x(n)) = 0.

y(n+1) = y(n−1)

import numpy as np
import matplotlib.pyplot as plt
plt.style.use("seaborn-poster")

def get_a_b(n):
    h = (np.pi/2-0) / n
    x = np.linspace(0, np.pi/2, n+1)
    # Get A
    A = np.zeros((n+1, n+1))
    A[0, 0] = 1
    A[n, n] = -2+4*h**2
    A[n, n-1] = 2
    for i in range(1, n):
        A[i, i-1] = 1
        A[i, i] = -2+4*h**2
        A[i, i+1] = 1
    # Get b
    b = np.zeros(n+1)
    for i in range(1, n+1):
        b[i] = 4*h**2*x[i]
    return x, A, b

x = np.pi/2
v = x - np.sin(2*x)
n_s = []
errors = []
for n in range(3, 100, 5):
    x, A, b = get_a_b(n)
    y = np.linalg.solve(A, b)
    n_s.append(n)
    e = v - y[-1]
    errors.append(e)
plt.figure(figsize = (10,8))
plt.plot(n_s, errors)
plt.yscale("log")
plt.xlabel("n gird points")
plt.ylabel("errors at x = $\pi/2$")

plt.savefig('ex226.png', dpi=72)
plt.show()

Last Updated on 2021-12-04 by gantovnik