Numerical Analysis I, Fall 2011 (http://www.math.nthu.edu.tw/˜ wangwc/)

Programming Homework Assignment

Assigned Dec 31, 2011. Due Jan 13, 2012.

1. The objective of this homework is to solve numerically the solution of the heat equation

u_{t}(x, y, t) = ν(u_{xx}(x, y, t) + u_{yy}(x, y, t)), (x, y) ∈ (0, 1)^{2}, t > 0  (1)

with initial and boundary conditions

u(x, y, 0) = sin(πx) sin(πy), u(x, 0, t) = u(x, 1, t) = u(0, y, t) = u(1, y, t) = 0, (2)

in order to get the solution u(x, y, t) at t = 1. Here ν is a material dependent parameter. We take ν = 0.1 in this problem.

The discrete approximation of u(x, y, t) is defined on (x_{i}, y_{j}, t^{n}) and is denoted by u_{i,j}^{n}, where x_{i} = iΔx, y_{j} = jΔy, t^{n} = nΔt. For simplicity, we take Δx = Δy = Δt = h = \frac{1}{N}. The numerical solution is obtained through an updating process, from {u_{i,j}^{n}}_{i,j=1,···N−1} to {u_{i,j}^{n+1}}_{i,j=1,···N−1} until one reaches n = N, t^{n} = t^{N} = 1, using the following formula

\begin{aligned}\frac{u_{i,j}^{n+1}−u_{i,j}^{n}}{Δt} = ν · \frac{1}{2}( & \frac{u_{i+1,j}^{n+1}−2u_{i,j}^{n+1}+u_{i−1,j}^{n+1}}{(Δx)^{2}} + \frac{u_{i,j+1}^{n+1}−2u_{i,j}^{n+1}+u_{i,j−1}^{n+1}}{(Δy)^{2}} \\ & \frac{u_{i+1,j}^{n}−2u_{i,j}^{n}+u_{i−1,j}^{n}}{(Δx)^{2}} + \frac{u_{i,j+1}^{n}−2u_{i,j}^{n}+u_{i,j−1}^{n}}{(Δy)^{2}} ) \end{aligned}  (3)

One can verify easily that the exact solution for (1), (2) is u^{e}(x, y, t) = e^{−2π^{2}νt} sin(πx) sin(πy). If all things are done correctly, you should see the error at t = 1 to scale like

max_{i,j=1,···N−1} |u^{e}(x_{i}, y_{j}, 1) − u_{i,j}^{N}| = O(Δx^{2}) + O(Δt^{2}) (= O(N^{−2}) in this problem) (4)

when you change N from 100 to 200, or from 200 to 400, etc.

(a) Start with the 1D case, where u^{e}(x, t) = e^{−π^{2}νt} sin(πx). Analyze possible approaches (ie. direct methods, iterative methods, etc) for solving the corresponding linear system and determine the quickest method and estimate the number of multiplication needed in terms of N. If you cannot do the 2D case, analyze, implement and hand in the 1D case instead for (minimal) partial credit.

(b) Continue on the 2D case. Analyze possible methods you learned from this class and explain why you chose your method. Hand in your analysis and your test result for (4) and your code (An email account for handing in the code will be announced later).

(c) (This one is a little harder and carries extra credit. Do it if time permits and be sure to consult me before and during the process) Solve numerically for

u_{xx}(x, y) + u_{yy}(x, y) − u^{3}(x, y) = sin(πx) sin(πy), (x, y) ∈ (0, 1)^{2}.

with the boundary condition

u(x, 0) = u(x, 1) = u(0, y) = u(1, y) = 0.

Now you have a nonlinear system to solve using, for example, the fixed point iteration or Newton’s method. A good initial guess for the fixed point/Newton iteration is given by u^{0}(x, y) = \frac{−1}{2π^{2}} sin(πx) sin(πy) (which is obtained by dropping the u^{3} term). In addition, there is another direct/iterative method needed for the linear system in either of the iteration. Analyze possible combination of methods mentioned above and explain the reason for your choices. Check that your answer is second order accurate (that is, the error is O(N^{−2})). Note that the exact solution u^{e}(x, y) is not known and you have to think about how to do this.

Remark:

Other possible discretization for (1) include

\frac{u_{i,j}^{n+1} − u_{i,j}^{n}}{Δt} = ν (\frac{u_{i+1,j}^{n+1} − 2u_{i,j}^{n+1} + u_{i−1,j}^{n+1}}{(Δx)^{2}} + \frac{u_{i,j+1}^{n+1} − 2u_{i,j}^{n+1} + u_{i,j−1}^{n+1}}{(Δy)^{2}} )  (5)

or

\frac{u_{i,j}^{n+1} − u_{i,j}^{n}}{Δt} = ν (\frac{u_{i+1,j}^{n} − 2u_{i,j}^{n} + u_{i−1,j}^{n}}{(Δx)^{2}} + \frac{u_{i,j+1}^{n} − 2u_{i,j}^{n} + u_{i,j−1}^{n}}{(Δy)^{2}} )  (6)

Both discretization give rise to error of O(Δx^{2}) + O(Δt) (= O(N^{−1}) for this problem). In addition, for (6) to be applicable, it is necessary that Δt ≤ CΔx^{2} (otherwise, you will get ’NaN’, which means ’not a number’, an indication for overflow, almost immediately).

Ideally, the best solution would be (3) together with a fast solver for the corresponding linear system. Such a faster solver (O(N^{2}) operations for the (N − 1)^{2} × (N − 1)^{2} matrix) does exist, but is beyond the scope of this course.

The next best method is probably (6) with Δt = CΔx^{2}. Since the purpose of this homework is to review your techniques for solving linear systems, we will not consider this option here.