{{ message }}

Instantly share code, notes, and snippets.

# ilciavo/BoundaryConditions.py

Last active Jun 4, 2020
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 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