Method of Manufactured Solutions for FDM in 1D

FDM_dirichlet.py

from __future__ import division
import numpy as np
from numpy import cos, sin, exp, zeros
from numpy.linalg import solve, norm
import matplotlib.pyplot as plt
from matplotlib import rcParams
 
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16 

def fdm_poisson(N, x, b_fun):    
    dx = x[1] - x[0]
    A = np.diag(-2*np.ones(N), 0)/dx**2 + np.diag(np.ones(N - 1), -1)/dx**2 \
    + np.diag(np.ones(N - 1), 1)/dx**2
    b = b_fun(x)
    u = zeros(N + 2)
    u[1:-1] = solve(A, b[1:-1])
    
    return u


def u_1(x):
    return (sin(x)**2 + x**2*exp(x))*(x - L/2)*(L/2 + x)


def f_1(x):    
    b = (-2*sin(x)**2 + 2*cos(x)**2 + x**2*exp(x) + 4*x*exp(x) 
        + 2*exp(x))*(x - L/2)*(L/2 + x) + 2*(2*cos(x)*sin(x) 
        + x**2*exp(x) + 2*x*exp(x))*(L/2 + x) + 2*(2*cos(x)*sin(x) 
        + x**2*exp(x) + 2*x*exp(x))*(x - L/2) + 2*sin(x)**2 + 2*x**2*exp(x)
    return b


def u_2(x):
    return (-sin(5*x)**2 - exp(-x**2) - x**2*exp(x))*(x - L/2.0)*(L/2.0 + x)


def f_2(x):
    b = (50*sin(5*x)**2 - 50*cos(5*x)**2 - 4*x**2*exp(-x**2)+2*exp(-x**2) 
        -x**2*exp(x) - 4*x*exp(x) - 2*exp(x))*(x - L/2.0)*(L/2.0 + x) \
        + 2*(-10*cos(5*x)*sin(5*x) + 2*x*exp(-x**2) - x**2*exp(x) 
        -2*x*exp(x))*(L/2.0 + x) + 2*(-10*cos(5*x)*sin(5*x) + 2*x*exp(-x**2)
        -x**2*exp(x) - 2*x*exp(x))*(x - L/2.0) - 2*sin(5*x)**2 - 2*exp(-x**2)\
        -2*x**2*exp(x)
    return b


L = 2
N_vec = np.array([5, 10, 50, 100, 500, 1000])
rel_error = np.zeros(N_vec.shape, dtype=float)
for cont, N in enumerate(N_vec):
    x = np.linspace(-0.5*L, 0.5*L, N + 2)
    u_fdm = fdm_poisson(N, x, f_2)
    u_anal = u_2(x)
    rel_error[cont] = norm(u_fdm - u_anal)/norm(u_anal)
    plt.plot(x, u_fdm, label='N=%i'%N)
    

plt.plot(x, u_anal, label='Analytic')
plt.xlabel(r"$x$")
plt.ylabel(r"$u(x)$")
plt.figure()
plt.loglog(N_vec, rel_error)
plt.xlabel(r"$N$")
plt.ylabel(r"Relative error")
plt.show()

FDM_periodic.py

from __future__ import division
import numpy as np
from numpy import cos, sin, zeros, pi
from numpy.linalg import solve, norm
import matplotlib.pyplot as plt
from matplotlib import rcParams
 
rcParams['font.family'] = 'serif'
rcParams['font.size'] = 16 

def fdm_poisson(N, x, b_fun):    
    dx = x[1] - x[0]
    A = np.diag(-2*np.ones(N), 0)/dx**2 + np.diag(np.ones(N-1), -1)/dx**2 \
    + np.diag(np.ones(N-1), 1)/dx**2
    A[0,N-1] = A[0,N-1] + 1
    A[N-1,0] = A[N-1,0] + 1
    b = b_fun(x)
    u = zeros(N)
    u[1:] = solve(A[1:,1:], b[1:])
    u[0] = u[-1]
    
    return u


def u_1(x):
    return 5*cos(5*pi*x)**2+sin(pi*x)


def f_1(x):    
    b = 250*pi**2*sin(5*pi*x)**2-250*pi**2*cos(5*pi*x)**2-pi**2*sin(pi*x)
    return b


L = 2
N_vec = np.array([50, 100, 500])
rel_error = np.zeros(N_vec.shape, dtype=float)
for cont, N in enumerate(N_vec):
    x = np.linspace(-0.5*L, 0.5*L, N)
    u_fdm = fdm_poisson(N, x, f_1)
    u_anal = u_1(x)
    rel_error[cont] = norm(u_fdm - u_anal)/norm(u_anal)
    plt.plot(x, u_fdm, label='N=%i'%N)
    

plt.plot(x, u_anal, 'k--', label='Analytic')
plt.xlabel(r"$x$")
plt.ylabel(r"$u(x)$")
plt.legend()
plt.figure()
plt.loglog(N_vec, rel_error)
plt.xlabel(r"$N$")
plt.ylabel(r"Relative error")
plt.show()

I thin poison is not well-posed problem with boundary conditions. You could try return (sin(x)2 + x2_exp(x))_(x - L/2)*(L/2 + x)-5