课程名称︰计算数学导论
课程性质︰数学系大三必修课
课程教师︰王伟仲
开课学院：理学院
开课系所︰数学系
考试日期（年月日）︰2015/1/14
考试时限（分钟）：3:30 ~ 6:00
试题 :
1. [25 points (3+3+6+10+3)] Consider the boundary value problem u'' + (x + 1)u'
2 -x -x
= 2u + (1-x )e , u(0) = -1, u(1) = 0, with exact solution u(x) = (x-1)e .
(a) If we want to solve the problem by using second-order finite difference,
how should we approximate the derivatives u''(x) and u'(x)? Just write down
your formula.
(b) Verify that the approximations proposed in (a) is truly of second-order
accuracy.
(c) If the interval [0,1] is divided into N regular subintervals and the
boundary value problem is discretized into Au = b where u = (u , u , u ,
0 1 2
..., u ), where u = -1 and u = 0. Form the matrix A and the vector b.
N 0 N
(d) [CODY] Write a program to solve the boundary value problem above. Input: h,
the length of the regular subintervals with h = 1/N. Output: e(u) = ∥ exact
sol. - numerical sol. ∥
2
(e) Complete the error table below and plot the error e(h) against h on a
log-log scale.
————————————————————————————————————
∣ h ∣ e(h) ∣ Rate of Convergence ∣
∣————— ∣———————————∣————————————————∣
∣ -1 ∣ ∣ ∣
∣ 10 ∣ ∣ ∣
∣————— ∣———————————∣————————————————∣
∣ -2 ∣ ∣ ∣
∣ 10 ∣ ∣ ∣
∣————— ∣———————————∣————————————————∣
∣ -3 ∣ ∣ ∣
∣ 10 ∣ ∣ ∣
∣————— ∣———————————∣————————————————∣
∣ -4 ∣ ∣ ∣
∣ 10 ∣ ∣ ∣
＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿＿
2. [10 points] [CODY] Write a program to implement the following method to
estimate the area of a quater circle with radius r = 0.5. That is , to estimate
r 2 2 r
I = ∫ √(r - x ) dx =∫ f(x) dx. Testing whether randomly generated points in
0 0
the square [0,r] ×[0,r] are inside or outside the quarter circle Ω(points on
the boundary are considered "inside") and multiplying the fraction inside by
2
r , the area of the rectangle.
Input: One N ×2 random array containing N points of random numbers in the
square
Output: Area of the quarter circle calculated
3. [15 points (5+5+5)]
(a) Determine the weights ω and points x of three-point Gaussian quadrature
i i
where i = 1,2,3.
b
(b) Suppose we apply the Composite Simpson's Rule to compute I(f) =∫ f(x) dx
a
on 2m regular subintervals where m is a positive integer. Given the error term
5
(x - x )
2j 2j-2 (4)
of Simpson's Rule in the interval [x , x ] is -————————f (ξ),
2j-2 2j 2880 j
where x < ξ < x and j = 1, 2, ..., m.
2j-2 j 2j
4
(b - a) h (4)
Prove the error term of Composite Simpson's Rule is -—————— f (ξ)
180
where a < ξ < b.
13.44 7 5 9
(c) If we want to compute I = ∫ (x - 2x + x - 1) dx with the accuracy of
-2.56
-8
10 , how many points in the interval [-2.56, 13.44] should we use at least in
order to achieve that accuracy in Gaussian Quadrature Rule and Composite
Simpson's Rule respectively?
4. [15 points (8+7)]
THE PCG ALGORITHM
[x,r] = pcg(A,M,b,x,tol)
Given symmetric positive defintie matrices A and M, a vector b, an initial
guess x, and a tolerance tol, compute an approximate solution to Ax = b.
T
Let r = b - Ax, solve Mz = r for z, and let γ = r z, p = z, and ρ = ∥r∥.
for k = 0, 1, ..., until ∥r∥/ρ < tol,
T
α=γ/(p Ap)
x = x + αp
r = r - αAp
Solve Mz = r for z.
^ T
γ = r z
^ ^
β = γ/γ, γ = γ
p = z + βp
end (for k)
(a) Go to "http://goo.gl/SA0EJH" to download the unfinished preconditioned
conjugate gradient (PCG) code and the test matrix A and the right hand side b.
Complete the code to solve the linear systems Ax = b by the following two ways:
(1) without using any preconditioner and (2) by using the diagonal
(0) T
preconditioner (i.e. M = diag(a)). Let the initial vector x = [0,0,...,0]
-8
and the tolerance tol = 10 . Plot two figures showing the number of iteration
versus the residual with or without preconditioning and upload them to ceiba
(save as .fig file).
(b) Modify the PCG code to implement the symmetric successive overrelaxation
preconditioning scheme (SSOR). That is, the preconditioner
ω 1 -1 1 T T
M = ——— (— D + L) D (— D + L) , where A = D + L + L , D is the diagonal
2-ω ω ω
of A, L is the strict lower triangular part of A. Your code should have the
inputs and outputs parameters [x,r] = pcg_ssor(A,b,omega,x,tol). You can use
MATLAB commands such as "\". Write down your code on the paper.
5. [15 points (9+3+3)]
(a) [CODY] Implement Euler's method to solve y'(t) = 3 - 2t - (1/2)y(t) for
t = [0,10] and y(0) = 1. The input is h and the output is t and y(t) such that
t is in [0,10].
(b) Use the code in (a) to solve the problem with h = 0.5. Plot the figure of t
versus computed y(t) and upload it to Ceiba (save as .fig file).
(c) Explain how you can use MATLAB ode23 to solve the initial value problem
defined in (a). Write down your code on paper.
6. [18 points (6+6+6)] Suppose the n ×n matrix A has eigenvalues λ , λ ,...
1 2
,λ ordered by |λ | > |λ | > |λ | ≧...≧ |λ |, with corresponding
n 1 2 3 n
(1) (2) (n)
eigenvectors v , v , ..., v that are linearly independent.
(a) Show that the Power method will converge to λ if we choose the initial
2
(0) (2) (3) (n)
vector x = β v + β v + ...+ β v .
2 3 n
(0)
(b) Show tha the initial vector x = (A -λ I)x satisfies the property given
1
(1) (2) (n)
in part (a). (Hint: Assume x = β v + β v + ...+ β v is a nonzero
1 2 n
vector.)
(c) Suppose you have known the first eigenpair, use (a) and (b) to design a
scheme to approximate λ .
2
7. [12 points] Categorize 6 numerical methods we have learned in this semester
into “Approximation”, “Iteration”, or “Equivalent Transformation”. In
each category, you should include 1 to 3 numerical methods. Explain why you
think your categorizations are meaningful.
8. [10 points] Use XMind (https://www.xmind.net/ ) to illustrate a concept map
of the topics we have learnt in this semester. You can pick the topic(s) you
are especially interested in.
=== End of problems. Good luck and have a happy winter break! ===