Final Project – Math 104B, Spring 20201 Instructor: Carlos J. García Cervera
1. (30 points) Consider the following partial differential equation for u on the square D = [0,1]2:
In (1),
−∆u = f(x,y), 0≤x,y≤1,
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)
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 Data.csv contains the experimental results. You can read the data using
import csv
T = []
P = []
with open(‘Data.csv’, mode=’r’, newline=”) as in_file:
in_reader = csv.reader(in_file , delimiter=’,’) for row in in_reader:
T = np.append(T,float(row[0])) P = np.append(P,float(row[1]))
(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)
(xi, yj−1)
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.