Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Poisson and Heat Equation with boundary conditions
#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()
@mgnisia

This comment has been minimized.

Copy link

@mgnisia mgnisia commented Jun 4, 2020

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

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment