Final Project – Math 104B, Spring 20201 Due on Friday, June 12th, 2020 Instructor: Carlos J. García Cervera
1. (30 points) Consider the following partial di erential equation for u on the square D = [0,1]2:
In (1),
u = f(x,y), 0x,y1,
u = 0 on@D. (1)
u = @2u + @2u, (2) @x2 @y2
and @D is the boundary of D:
@D = {0, 1} ⇥ [0, 1] [ [0, 1] ⇥ {0, 1}. (3)
We construct a two-dimensional grid: Let h = 1/n, and set xi = ih, i=0,1,…,n,
yj = jh, j=0,1,…,n. (4) We discretize the equation in a similar way to what we did in lecture:
@2u(xi,yj) = u(xi +h,yj) 2u(xi,yj)+u(xi h,yj) +O(h2), (5) @x2 h2
and
@2u(xi,yj) = u(xi,yj +h) 2u(xi,yj)+u(xi,yj h) +O(h2). (6)
@y2 h2
1All course materials (class lectures and discussions, handouts, homework assignments, examinations, web materials, etc) and the intellectual content of the course itself are protected by United States Federal Copyright Law, and the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lecture notes and all other course materials without the prior written permission of the instructor.
1
2
Therefore we set up the system
ui+1,j 2ui,j+ui 1,j ui,j+1 2ui,j+ui,j 1 =fi,j, i,j=1,2,…,n. h2 h2 (7)
(a). (5 points) Construct the function f (x, y) so that u(x, y) = sin(16⇡x(1 x)y(1 y)) solves the equation (see Figure 1).
Figure 1: Solution to equation (1).
(b). (25 points) Solve system (7) using the Conjugate Gradient Method. Set the tolerance to 10 8, and use n = 10, 20, 40, 80, 160, 320. Plot the timings and the number of iterations, and determine the scaling of each in terms of n. Verify numerically that the approximation is second order accurate.
Implementaion Tips:
(a). for a simple implementation you can define the solution u (and all arrays in the Conjugate Gradient Method) as two-dimensional arrays, just like it is defined in (7).
3
(xi, yj+1)
(xi,yj)
(xi+1,yj)
(xi 1,yj)
(xi, yj 1)
Figure 2: Indexing of the nodes in the two-dimensional grid.
(b). Do not construct the full matrix associated to system (7). Instead, there are other things you can do, for example you can define the matrix using sparse linear algebra, or simply define a function to evaluate the matrix/vector product directly:
(Ap)ij = pi+1,j 2pi,j + pi 1,j pi,j+1 2pi,j + pi,j 1 . (8) h2 h2
2. (30 points) A physical quantity P is known to depend on the temperature T . The file contains the experimental results. You can read the data using
(a). (20 points) Use the Least Squares method with the Householder QR to fit the data by a polynomial model of the form
P(T)=c0 +c1T +c2T2 +…+ckTk, (9)
4
for k = 1, 2, 3, 4, 5. Compute the modeling error in each case. Plot the values
(Ti, Pi) and the results from each model for comparison.
(b). (10 points) Do a log-log plot of the same values. Explain why the log-log plot
indicates that P (T ) = aT b for some a and b. Use the Least Squares method and solve the Normal Equations to estimate the values of a and b, and do a regular plot and a log-log plot with the values (Ti, Pi) and the values (Ti, aTib) to show how good your approximation is. Compute the modeling error as well.
Extra Credit – Math 104B, Spring 20201 Due on Friday, June 12th, 2020 Instructor: Carlos J. García Cervera
1. (40 points) When discretizing the di erential equation u00(x) = f(x), x 2 [0, 1]
u(0) = u(1) = 0, (1) using second order centered di erences, we obtained the tridiagonal matrix
0B2 1 0 0 … 0 01C B 1 2 1 0 … 0 0 C B0 1 2 1 … 0 0C
A= 1 B . . … … … . . C. (2)
h2 B 0 0 0 1 2 1 0 C B@0 0 0 0 1 2 1CA
0 0 0 0 0 1 2 In components, remember that this means
ui+1 2ui+ui 1 =f(xi) i=1,2,…,n, (3) h2
where h = 1/(n + 1).
(a). (10 points) For each k = 1,2,…,n, show that the vector u(k) given by
u(k) = sin✓ ⇡ki ◆, i = 1,2,…,n (4) i n+1
1All course materials (class lectures and discussions, handouts, homework assignments, examinations, web materials, etc) and the intellectual content of the course itself are protected by United States Federal Copyright Law, and the California Civil Code. The UC Policy 102.23 expressly prohibits students (and all other persons) from recording lectures or discussions and from distributing or selling lecture notes and all other course materials without the prior written permission of the instructor.
1
2
is an eigenvector of the matrix A, and determine the corresponding eigenvalue
k.
(b). (5 points) The Jacobi iteration matrix is defined as
TJ =D 1(L+U), (5)
where A = L + D + U. Show that the vectors (12) are also eigenvectors of TJ and find its eigenvalues. Hint: Prove that since D = dI (constant multiple of the identity),
TJ x = μx () Ax = d(μ + 1)x. (6)
(c). (5 points) Determine the spectral radius of TJ,
⇢(TJ) = max |μk|. (7) k=1,…,n
.
(d). (5 points) The Jacobi iteration can be written as
x(k+1) = TJx(k) + c. (8) From the previous steps, we know that TJ is symmetric and diagonalizable.
Use this fact to show that if x⇤ is the (unique) fixed point of (8), then
kx(k) x⇤k2 ⇢(TJ)kkx(0) x⇤k2 (9)
(e). (5 points) Use formula (9) and the spectral radius obtained earlier to show that the method converges, and to estimate the number of iterations necessary for the error to be less than a given ✏ as a function of the number of grid points used, n. You should end up with a formula of the form Iter(n) = O(n↵) for some ↵. Hint: At the end you might need to use Taylor to estimate ↵.
(f). (5 points) To verify the value of ↵ numerically, fix ✏ = 10 6. Consider an initial random vector x(0), and solve the system of equations
Ax = 0 (10)
3
using Jacobi’s method. Use the values n = 10, 20, 40, 80, 160, 320. Do a log-log plot of the number of Jacobi iterations necessary for the error to satisfy
kx(k)k2 ✏kx(0)k2. (11)
How does the computed number of iterations compare with the theoretical one
obtained earlier?
(g). (5 points) Repeat the previous part with Gauss-Seidel’s method. How much
faster is it?
2. (20 points) In this problem, we analyze the Gauss-Seidel method in the same way
we studied Jacobi’s method above:
(a). (10 points) For each k = 1,2,…,n, show that the vector u(k) given by
(k) ✓ ✓k⇡◆◆i ✓⇡ki◆
ui = cos n+1 sin n+1 , i=1,2,…,n (12)
is an eigenvector of the corresponding Gauss-Seidel iteration matrix, and deter- mine the corresponding eigenvalue k . Hint: Verify directly that (L+D)u(k) = kUu(k) for some k.
(b). (5 points) Use the previous part to show that ⇢(TGS) = ⇢(TJ )2, where TGS and TJ are the iteration matrices for the Gauss-Seidel and the Jacobi methods, respectively.
(c). (5 points) Explain why if one uses the Gauss-Seidel method for this equation, one only needs half the number of iterations of the Jacobi method to reach a tolerance ✏.