Poisson and Heat Equation with boundary conditions

BoundaryConditions.py

#poisson solver
#u_xx = b(x)
%matplotlib inline 
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'] = 15
 
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 = solve(A,b)
    u[1:] = solve(A[1:,1:], b[1:])
    u[0] = u[-1]
    return u

#Analytical solution
def u_1(x):
    return 5*cos(5*pi*x)**2+sin(pi*x)-5
 
#right hand side
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)
    b = -250*pi**2*cos(10*pi*x)-pi**2*sin(pi*x)
    return b
 
 
L = 2
N_vec = np.array([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)
    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()

from numpy import linspace, exp, diag, ones, append,delete,dot
# Solve: u_t = u_xx + b
# with periodic boundary conditions

# analytical solution:
def uRef(t,x): 
    return exp(-t)*cos(5*pi*x)
    
def b(t,x):
    return (25*pi**2-1)*exp(-t)*cos(5*pi*x)

# grid
N = 30
Tend = 2.
x = linspace(-1,1,N+1)

# leave off 1 point so initial condition is periodic
# (doesn't have a duplicate point)
x = delete(x,-1)
uWithMatrix = uRef(0,x)
uNoMatrix = uRef(0,x)
uImplicit = uRef(0,x)

dx = x[1]-x[0]
dt = dx**2/2

#Laplacian matrix:
A = diag(-2*ones(N), 0)+ diag(ones(N-1), -1)+diag(np.ones(N-1), 1)
A[0,N-1] = 1
A[N-1,0] = 1
A = A/dx**2

iLeft = append(len(x)-1, range(0,len(x)-1))
iCenter = range(0,len(x))
iRight = append(range(1,len(x)), 0)

for t in linspace(0,Tend-dt,1/dt):
    b_val = b(t,x)
    uImplicit = solve(1-A*dt, dt*b_val)
    uWithMatrix = uWithMatrix + dt*(dot(A,uWithMatrix) + b(t,x))
    uNoMatrix = uNoMatrix + dt*((uNoMatrix[iLeft] - 2*uNoMatrix[iCenter] + uNoMatrix[iRight])/dx**2 + b(t,x))
    
uAnal = uRef(Tend-dt,x)

fig = plt.figure(figsize=(15,10))
plt.plot(x, uAnal, 'k--', label='Analytical Solution')
plt.plot(x, uWithMatrix, 'b-o', label='Matrix Solution')
plt.plot(x,uImplicit, 'r-o', label='Implicit')
plt.plot(x, uNoMatrix, 'g-o', label='No Matrix Solution')
plt.xlabel(r"$x$")
plt.ylabel(r"$u(x)$")
plt.legend()

Thank you for this example! :) One slight modification was necessary to make it work again: replace in l.94

linspace(0,Tend-dt,1/dt)

with

linspace(0,Tend-dt,int(1/dt))

otherwise numpy throws a TypeError as the number of elements is not an integer