Num Methods: Finite Differences in 1D

30 days ago by solin

Finite Differences for 1D 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 1D with Dirichlet and Neumann boundary conditions, and how to solve a time-dependent problem du/dt - ddu/dxx = f.

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

Equation -u'' = f with Dirichlet boundary conditions (solution values prescribed) at both interval endpoints.
# This function solves the equation -u'' = f in interval (a, b) with # Dirichlet conditions on both ends. The interval is subdivided into # 'n' subintervals. def fdm_solve_1d_1(a, b, n, bc_dir_left, bc_dir_right, f): # matrix size size = n-1 # equidistant grid h = (b-a)/float(n) x_array = arange(a+h, b-h/2., h) # right-hand side rhs = [f(xx)*h*h for xx in x_array] # incorporate boundary conditions rhs[0] += bc_dir_left rhs[-1] += bc_dir_right # assemble matrix print "Assembling." I, J, V = generate_tridiagonal_coo_matrix(size, 0, 2., -1.) # solve the matrix problem print "Solving." m = sparse.coo_matrix((V,(I,J)),shape=(size,size)) m = m.tocsr() rhs = array(rhs) #sol, res = gmres(m, rhs, tol=1.e-8) # also possible: cg, cgs, qmr, gmres, bicg, bicgstab, ... sol = spsolve(m, rhs) # plot solution x_array_ext = array([a] + list(x_array) + [b]) sol_ext = array([bc_dir_left] + list(sol) + [bc_dir_right]) print "Plotting." clf() axis('equal') label = "solution" plot(x_array_ext, sol_ext, "g-", label=label) legend() grid(True) savefig("a.png") # Define load function def f(x): #return -1. if x > (a + b)/2.: return 1. else: return -1. # Problem parameters a = 0. b = 1. n = 10 bc_dir_left = 1. bc_dir_right = 1 # Solve by FDM fdm_solve_1d_1(a, b, n, bc_dir_left, bc_dir_right, f) 
       

Case 2

Equation -u'' = f with a Dirichlet boundary condition (solution value) on the left and Neumann boundary condition (solution derivative) on the right.
# This function solves the equation -u'' = f in interval (a, b) with # a Dirichlet condition on the left and a Neumann condition on the right. # The interval is subdivided into 'n' subintervals. def fdm_solve_1d_2(a, b, n, bc_dir_left, bc_neum_right, f): # matrix size size = n # equidistant grid h = (b-a)/float(n) x_array = arange(a+h, b+h/2., h) # right-hand side rhs = [f(xx)*h*h for xx in x_array] # assemble matrix print "Assembling." I, J, V = generate_tridiagonal_coo_matrix(size, 0, 2., -1.) # incorporate boundary conditions # Dirichlet left: rhs[0] += bc_dir_left # Neumann right: #m[n-1][n-2] = -1. I.append(size-1) J.append(size-2) V.append(1. - 1.) #m[n-1][n-1] = 1. I.append(size-1) J.append(size-1) V.append(-2. + 1.) rhs[-1] += h*bc_neum_right # solve the matrix problem print "Solving." m = sparse.coo_matrix((V,(I,J)),shape=(size,size)) m = m.tocsr() rhs = array(rhs) #sol, res = gmres(m, rhs, tol=1.e-8) # also possible: cg, cgs, qmr, gmres, bicg, bicgstab, ... sol = spsolve(m, rhs) # plot solution x_array_ext = array([a] + list(x_array)) sol_ext = array([bc_dir_left] + list(sol)) print "Plotting." clf() axis('equal') label = "solution" plot(x_array_ext, sol_ext, "g-", label=label) legend() grid(True) savefig("a.png") # Define load function def f(x): return -1. # Problem parameters a = 0. b = 1. n = 100 bc_dir_left = 1. bc_neum_right = 0 # Solve by FDM fdm_solve_1d_2(a, b, n, bc_dir_left, bc_neum_right, f) 
       

Case 3

Equation -u'' = f with a Dirichlet boundary condition (solution value) on the left and Newton boundary condition (combination of solution value and derivative) on the right.
# This function solves the equation -u'' = f in interval (a, b) with # a Dirichlet condition on the left and a Newton condition c1*u(x) + c2*du/dx = c3 # on the right. The interval is subdivided into 'n' subintervals. def fdm_solve_1d_3(a, b, n, bc_dir_left, bc_newton_c1, bc_newton_c2, bc_newton_c3, f): # matrix size size = n # equidistant grid h = (b-a)/float(n) x_array = arange(a+h, b+h/2., h) # right-hand side rhs = [f(xx)*h*h for xx in x_array] # assemble matrix print "Assembling." I, J, V = generate_tridiagonal_coo_matrix(size, 0, 2., -1.) # incorporate boundary conditions # Dirichlet left: rhs[0] += bc_dir_left # Newton right: # m[n-1][n-2] = -bc_newton_c2 I.append(size-1) J.append(size-2) V.append(1. - bc_newton_c2) #m[n-1][n-1] = h*bc_newton_c1 + bc_newton_c2 I.append(size-1) J.append(size-1) V.append(-2. + h*bc_newton_c1 + bc_newton_c2) rhs[-1] = h*bc_newton_c3 # solve the matrix problem print "Solving." m = sparse.coo_matrix((V,(I,J)),shape=(size,size)) m = m.tocsr() rhs = array(rhs) #sol, res = gmres(m, rhs, tol=1.e-8) # also possible: cg, cgs, qmr, gmres, bicg, bicgstab, ... sol = spsolve(m, rhs) # plot solution x_array_ext = array([a] + list(x_array)) sol_ext = array([bc_dir_left] + list(sol)) print "Plotting." clf() axis('equal') label = "solution" plot(x_array_ext, sol_ext, "g-", label=label) legend() grid(True) savefig("a.png") # Define load function def f(x): return -1. #if x > (a + b)/2.: return -1. #else: return 0 # Problem parameters a = 0. b = 1. n = 100 bc_dir_left = 1. bc_newton_c1 = 0. bc_newton_c2 = 1. bc_newton_c3 = 0. # Solve by FDM fdm_solve_1d_3(a, b, n, bc_dir_left, bc_newton_c1, bc_newton_c2, bc_newton_c3, f) 
       

Case 4

Time-dependent equation du/dt - ddu/dxx = f(x,t) with Dirichlet boundary conditions, discretized in time using implicit Euler method.
# This function solves the equation du/dt - ddu/dxx = f(x,t) in a spatial interval (a, b) # and temporal interval (0, T), with Dirichlet conditions on both ends and an initial # condition f_init(x). Note: the initial condition must be compatible with the boundary # conditions. The spatial interval is subdivided into 'nx' subintervals, the time interval # into 'nt' time steps. from numpy import pi def fdm_solve_1d_4(a, b, T, nx, nt, bc_dir_left, bc_dir_right, f, f_init): # matrix size size = nx-1 # equidistant grid in space h = (b-a)/float(nx) x_array = arange(a+h, b-h/2., h) # time step dt = T/float(nt) # time stepping loop current_time = 0. # grid for plotting purposes x_array_ext = array([a] + list(x_array) + [b]) # setting initial condition sol = [f_init(xx) for xx in x_array] sol_ext = array([bc_dir_left] + list(sol) + [bc_dir_right]) # plotting initial condition clf() axis('equal') label = "solution, t = " + str(current_time) plot(x_array_ext, sol_ext, "g-", label=label) legend() grid(True) filename = "a_" + str(current_time) + ".png" savefig(filename) # assemble matrix print "Assembling." I, J, V = generate_tridiagonal_coo_matrix(size, 0, 2. + h*h/dt, -1.) m = sparse.coo_matrix((V,(I,J)),shape=(size,size)) m = m.tocsr() # time stepping loop while(1): current_time += dt if current_time > T + dt/2.: return print "Time = ", float(current_time) # right-hand side rhs = [f(xx, current_time)*h*h for xx in x_array] for i in range(nx-1): rhs[i] += h*h/dt*sol[i] # incorporate boundary conditions rhs[0] += bc_dir_left rhs[-1] += bc_dir_right # solve the matrix problem print "Solving." rhs = array(rhs) #sol, res = gmres(m, rhs, tol=1.e-8) # also possible: cg, cgs, qmr, gmres, bicg, bicgstab, ... sol = spsolve(m, rhs) # plot solution sol_ext = array([bc_dir_left] + list(sol) + [bc_dir_right]) print "Plotting." clf() axis('equal') label = "solution, t = " + str(current_time) plot(x_array_ext, sol_ext, "g-", label=label) legend() grid(True) filename = "a_" + str(current_time) + ".png" savefig(filename) # Define load function def f(x, t): return sin(x)*(cos(t) + sin(t)); # Define initial condition def f_init(x): return 0. # Problem parameters a = 0. b = 2*pi T = pi nx = 100 nt = 25 bc_dir_left = 0. bc_dir_right = 0. # Solve by FDM fdm_solve_1d_4(a, b, T, nx, nt, bc_dir_left, bc_dir_right, f, f_init)