Num Methods: Finite Differences in 2D

30 days ago by solin

Finite Differences for 2D Problems

Finite Difference Method (FDM) is the most basic numerical method to solve partial differential equations. This worksheet shows you how to solve a Poisson equation -u'' = f in 2D with zero and nonzero Dirichlet boundary conditions.

The following active input text box contains the main code. You need to click on it to see an "Evaluate" link appear at its bottom. Click on this link to have the code sent to our computer and interpreted. Then you can evaluate/modify/evaluate the examples in the other text boxes in the same way. Have fun and let us know at femhub@googlegroups.com with any problems or questions.

The Code:
# Import symbolics, plotting, and linear algebra functions. from pylab import plot, savefig, grid, legend, clf, pcolor, spy, axis from numpy import arange, zeros, dot, array from numpy.linalg import norm from sympy import (var, pi, factorial, sin, cos, log, exp, tan, atan, asin, acos, sinh, cosh, tanh, asinh, acosh, lambdify) from scipy import shape, eye, tril, triu, diag, sparse from scipy.linalg import solve from scipy.sparse.linalg import cg, cgs, qmr, gmres, bicg, bicgstab from scipy.sparse.linalg.dsolve import spsolve from time import time # size...matrix size # o...offset of side-diagonal def generate_tridiagonal_full_matrix(size, o, a, b): m = zeros((size, size)) for i in range(0, size): for j in range(0, size): if i == j: m[i][i] = a if j == i + 1 + o: m[i][j] = b m[j][i] = b return m # size...matrix size # o...offset of side-diagonal def generate_tridiagonal_coo_matrix(size, o, a, b): I = [] J = [] V = [] for i in range(0, size): for j in range(0, size): if i == j: I.append(i) J.append(i) V.append(a) if j == i + 1 + o: I.append(i) J.append(j) V.append(b) I.append(j) J.append(i) V.append(b) return I, J, V 
       
 

Case 1

Two-dimensional example -Laplace u = f with nonzero Dirichlet boundary conditions in a square (0,a)x(0,b). The nonzero Dirichlet boundary conditions are given through a function g(x,y,a,b).
import numpy # Solve equation -Laplace u = f with Dirichlet boundary # conditions in the square (0, a) x (0, b). def fdm_solve_2d_1(a, b, nx, ny, f, g): hx = a/float(nx) hx2 = hx*hx hy = b/float(ny) hy2 = hy*hy size = (nx-1)*(ny-1) I = [] J = [] V = [] rhs = numpy.zeros(size) print "Assembling." for m1 in range(1, nx): for m2 in range(1, ny): if inside(m1, m2, nx, ny): # matrix entry i_c = (m1-1)*(ny-1) + (m2-1) I.append(i_c) J.append(i_c) V.append(2./hx2 + 2./hy2) x_phys = m1*hx y_phys = m2*hy coeff = -1./hx2 if away_from_left_bdy(m1, m2, nx, ny): i_l = i_c - (ny-1) I.append(i_c) J.append(i_l) V.append(coeff) else: rhs[i_c] = rhs[i_c] - coeff * g(x_phys-hx, y_phys, a, b) coeff = -1./hx2 if away_from_right_bdy(m1, m2, nx, ny): i_r = i_c + (ny-1) I.append(i_c) J.append(i_r) V.append(coeff) else: rhs[i_c] = rhs[i_c] - coeff * g(x_phys+hx, y_phys, a, b) coeff = -1./hy2 if away_from_bottom_bdy(m1, m2, nx, ny): i_b = i_c - 1 I.append(i_c) J.append(i_b) V.append(coeff) else: rhs[i_c] = rhs[i_c] - coeff * g(x_phys, y_phys-hy, a, b) coeff = -1./hy2 if away_from_top_bdy(m1, m2, nx, ny): i_t = i_c + 1 I.append(i_c) J.append(i_t) V.append(coeff) else: rhs[i_c] = rhs[i_c] - coeff * g(x_phys, y_phys+hy, a, b) # rhs entry rhs[i_c] = rhs[i_c] + f(x_phys, y_phys) # creating CSR matrix from arrays I, J, V m = sparse.coo_matrix((V,(I,J)),shape=(size,size)) m = m.tocsr() rhs = array(rhs) print "Solving." t = time() #sol, res = gmres(m, rhs, tol=1.e-8) # also possible: cg, cgs, qmr, gmres, bicg, bicgstab, ... sol = spsolve(m, rhs) print "Matrix solve took", time()-t, "s." # Plotting print "Plotting." from mpl_toolkits.mplot3d import Axes3D from pylab import figure, savefig import numpy as np from numpy import zeros, sin from matplotlib import cm fig = figure() ax = Axes3D(fig) X, Y = np.mgrid[0:(nx+1)*hx:hx, 0:(ny+1)*hy:hy] m = zeros((nx+1, ny+1), dtype="float") # defining values at Cartesian grid points for m1 in range(nx+1): for m2 in range(ny+1): if m1 > 0 and m2 > 0 and m1 < nx and m2 < ny: i_c = (m1-1)*(ny-1) + (m2-1) m[m1][m2] = sol[i_c] else: x_phys = m1*hx y_phys = m2*hy m[m1][m2] = g(x_phys, y_phys, a, b) # putting them into the array Z for plotting Z = m ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet) savefig("surface1.png") # Problem parameters # right-hand side def f(x,y): return 0. # Dirichlet boundary conditions def g(x,y,a,b): return x*y # domain edges and division a = 1. b = 2. nx = 10 ny = 20 # define the domain boundary def inside(m1, m2, nx, ny): return m1 >= 1 and m1 <= nx-1 and m2 >= 1 and m2 <= ny-1 def away_from_left_bdy(m1, m2, nx, ny): return m1 > 1 def away_from_right_bdy(m1, m2, nx, ny): return m1 < nx-1 def away_from_bottom_bdy(m1, m2, nx, ny): return m2 > 1 def away_from_top_bdy(m1, m2, nx, ny): return m2 < ny-1 # Solve by FDM fdm_solve_2d_1(a, b, nx, ny, f, g) 
       

Case 2

Two-dimensional example -Laplace u = f with (generally nonzero) Dirichlet boundary conditions in a general subdomain of the square (0,a)x(0,b). The Dirichlet boundary conditions are given through a function g(x,y,a,b).
import numpy # calculate the number of interior grid points def calc_matrix_size(nx, ny): size = 0 for m1 in range(1, nx): for m2 in range(1, ny): if inside(m1, m2, nx, ny): size = size + 1 return size # calculate the index of an interior grid point # (asserts that the grid point lies inside) def calc_index(m1_, m2_, nx, ny): i = 0 # interior grid point number (= equation number) for m1 in range(1, nx): for m2 in range(1, ny): if inside(m1, m2, nx, ny): if m1 == m1_ and m2 == m2_: return i i = i + 1 # Solve equation -Laplace u = f with Dirichlet boundary # conditions in a subdomain of the square (0, a) x (0, b). def fdm_solve_2d_2(a, b, nx, ny, f, g): hx = a/float(nx) hx2 = hx*hx hy = b/float(ny) hy2 = hy*hy size = calc_matrix_size(nx, ny) print "matrix size =", size I = [] J = [] V = [] rhs = numpy.zeros(size) print "Assembling." i_c = 0 # interior grid point number (= equation number) for m1 in range(1, nx): for m2 in range(1, ny): if inside(m1, m2, nx, ny): # matrix entry I.append(i_c) J.append(i_c) V.append(2./hx2 + 2./hy2) x_phys = m1*hx y_phys = m2*hy coeff = -1./hx2 if away_from_left_bdy(m1, m2, nx, ny): i_l = calc_index(m1-1, m2, nx, ny) I.append(i_c) J.append(i_l) V.append(coeff) else: rhs[i_c] = rhs[i_c] - coeff * g(x_phys-hx, y_phys, a, b) coeff = -1./hx2 if away_from_right_bdy(m1, m2, nx, ny): i_r = calc_index(m1+1, m2, nx, ny) I.append(i_c) J.append(i_r) V.append(coeff) else: rhs[i_c] = rhs[i_c] - coeff * g(x_phys+hx, y_phys, a, b) coeff = -1./hy2 if away_from_bottom_bdy(m1, m2, nx, ny): i_b = i_c - 1 I.append(i_c) J.append(i_b) V.append(coeff) else: rhs[i_c] = rhs[i_c] - coeff * g(x_phys, y_phys-hy, a, b) coeff = -1./hy2 if away_from_top_bdy(m1, m2, nx, ny): i_t = i_c + 1 I.append(i_c) J.append(i_t) V.append(coeff) else: rhs[i_c] = rhs[i_c] - coeff * g(x_phys, y_phys+hy, a, b) # rhs entry rhs[i_c] = rhs[i_c] + f(x_phys, y_phys) # increasing active grid point number i_c = i_c + 1 # creating CSR matrix from arrays I, J, V m = sparse.coo_matrix((V,(I,J)),shape=(size,size)) m = m.tocsr() rhs = array(rhs) print "Solving." t = time() #sol, res = gmres(m, rhs, tol=1.e-8) # also possible: cg, cgs, qmr, gmres, bicg, bicgstab, ... sol = spsolve(m, rhs) print "Matrix solve took", time()-t, "s." # Plotting print "Plotting." from mpl_toolkits.mplot3d import Axes3D from pylab import figure, savefig import numpy as np from numpy import zeros, sin from matplotlib import cm fig = figure() ax = Axes3D(fig) X, Y = np.mgrid[0:(nx+1)*hx:hx, 0:(ny+1)*hy:hy] m = zeros((nx+1, ny+1), dtype="float") # defining values at Cartesian grid points i_c = 0 for m1 in range(nx+1): for m2 in range(ny+1): if inside(m1, m2, nx, ny): m[m1][m2] = sol[i_c] i_c = i_c + 1 else: if on_boundary(m1, m2, nx, ny): x_phys = m1*hx y_phys = m2*hy m[m1][m2] = g(x_phys, y_phys, a, b) else: m[m1][m2] = 0 # putting them into the array Z for plotting Z = m ax.plot_surface(X, Y, Z, rstride=1, cstride=1, cmap=cm.jet) savefig("surface1.png") # Problem parameters # right-hand side def f(x,y): return 1. # Dirichlet boundary conditions def g(x,y,a,b): #return x*y return 0 # outer rectangle measures a = 1. b = 1. # divisions nx = 60 ny = 60 # deleted cells in the lower-left corner dx = 30 dy = 30 # define the domain boundary def on_boundary(m1, m2, nx, ny): if m1 == 0: return True if m1 == nx: return True if m2 == 0: return True if m2 == ny: return True if m1 == dx and m2 <= dy: return True if m1 <= dx and m2 == dy: return True return False def inside(m1, m2, nx, ny): if m1 < 1 or m1 > nx-1 or m2 < 1 or m2 > ny-1: return False if m1 <= dx and m2 <= dy: return False return True def away_from_left_bdy(m1, m2, nx, ny): if m1 > dx + 1: return True if m1 > 1 and m2 > dy: return True return False def away_from_right_bdy(m1, m2, nx, ny): return m1 < nx - 1 def away_from_bottom_bdy(m1, m2, nx, ny): if m2 > dy + 1: return True if m1 > dx and m2 > 1: return True return False def away_from_top_bdy(m1, m2, nx, ny): return m2 < ny-1 # Solve by FDM fdm_solve_2d_2(a, b, nx, ny, f, g)