師大
Tsung-Ming Huang
Department of Mathematics National Taiwan Normal University, Taiwan
March 11, 2009
師大
Outline
1 Model problem
2 Center difference discretization
師大
Consider the Dirichlet boundary-value problem:
−∆u ≡ −uxx− uyy = 2π2sin πx sin πy, for (x, y) ∈ Ω, (1) u(x, y) = 0 (x, y) ∈ ∂Ω,
for Ω := {x, y|0 < x, y < 1} ⊆ R2 with boundary ∂Ω, which has the exact solution
u(x, y) = sin πx sin πy.
師大
師大
To solve (1) by means of a difference methods, one replaces the differential operator by a difference operator. Let
Ωh := {(xi, yi)|i, j = 1, . . . , n},
∂Ωh := {(xi, 0), (xi, 1), (0, yj), (1, yj)|i, j = 0, 1, . . . , n + 1}, where xi = ih, yj = jh, i, j = 0, 1, . . . , n + 1, h := n+11 , n ≥ 1, is an integer.
師大
From the Taylor’s theorem, we have u(xi+ h) = u(xi) + u0(xi)h + h2
2 u00(xi) +h3
6 u000(xi) + h4
24u(4)(ξ1) u(xi− h) = u(xi) − u0(xi)h + h2
2 u00(xi) −h3
6 u000(xi) + h4
24u(4)(ξ2), where ξ1is between xi and xi+ hand ξ2 is between xiand
xi− h. Hence
u00(xi) = u(xi+ h) − 2u(xi) + u(xi− h)
h2 −h2
12u(4)(ξ)
= u(xi+1) − 2u(xi) + u(xi−1)
h2 −h2
12u(4)(ξ), where ξ is between xi− h and xi+ h.
師大
Similarly,
∂2u
∂x2(xi, yj) = u(xi+1, yj) − 2u(xi, yj) + u(xi−1, yj)
h2 −h2
12
∂4u
∂x4(ξi, yj),
∂2u
∂y2(xi, yj) = u(xi, yj+1) − 2u(xi, yj) + u(xi, yj−1)
h2 −h2
12
∂4u
∂x4(xi, ηj), where ξi∈ (xi−1, xi+1)and ηj ∈ (yj−1, yj+1). It implies that
∂2u
∂x2(xi, yj) + ∂2u
∂y2(xi, yj)
= u(xi, yj−1) + u(xi−1, yj) − 4u(xi, yj) + u(xi+1, yj) + u(xi, yj+1) h2
−h2 12
∂4u
∂x4(ξi, yj) + ∂4u
∂x4(xi, ηj)
.
師大
Let uij denote an approximated value of function u at the grid point (xi, yj)for i, j = 1, . . . , n + 1. Then
−uxx(xi, yj) − uyy(xi, yj)
≈ −ui,j−1− ui−1,j+ 4ui,j− ui+1,j− ui,j+1 h2
with an error O(h2)and the equation
−uxx(xi, yj) − uyy(xi, yj) = 2π2sin πxisin πyj ≡ fij can be replaced by the following equation
−ui,j−1− ui−1,j+ 4ui,j− ui+1,j− ui,j+1
h2 = fij (2)
for i, j = 1, . . . , n.
師大
For j = 1, we have
−u1,0− u0,1+ 4u1,1− u2,1− u1,2 = h2f1,1, (3a)
−u2,0− u1,1+ 4u2,1− u3,1− u2,2 = h2f2,1, (3b) ...
−un−1,0− un−2,1+ 4un−1,1− un,1− un−1,2 = h2fn−1,1,(3c)
−un,0− un−1,1+ 4un,1− un+1,1− un,2 = h2fn,1. (3d) By the boundary condition, it holds that
u1,0 = u2,0= · · · = un,0= 0, (4a)
u0,1 = un+1,1= 0. (4b)
師大
Substituting (4) into (3), we get
4u1,1− u2,1−u1,2 = h2f1,1, (5a)
−u1,1+ 4u2,1− u3,1−u2,2 = h2f2,1, (5b) ...
−un−2,1+ 4un−1,1− un,1−un−1,2 = h2fn−1,1, (5c)
−un−1,1+ 4un,1−un,2 = h2fn,1. (5d) Let, for j = 1, . . . , n,
u:,j =
u1,j
u2,j ... un,j
, f:,j =
f1,j
f2,j ... fn,j
, A1=
4 −1
−1 . .. ...
. .. ... −1
−1 4
∈ Rn×n.
師大
Then (5) can be rewritten as following matrix form:
A1 −In
u:,1
u:,2
= h2f:,1.
For j = 2, . . . , n − 1, using u0,j = un+1,j = 0, we have
−u1,j−1+ 4u1,j− u2,j−u1,j+1 = h2f1,j,
−u2,j−1− u1,j+ 4u2,j− u3,j−u2,j+1 = h2f2,j, ...
−un−1,j−1− un−2,j+ 4un−1,j− un,j−un−1,j+1 = h2fn−1,j,
−un,j−1− un−1,j+ 4un,j−un,j+1 = h2fn,j. Above equations can be represented as following matrix form:
−In A1 −In
u:,j−1
u:,j
u:,j+1
= h2f:,j.
師大
For j = n, using u1,n+1= u2,n+1= un,n+1 = 0, we have
−u1,n−1+ 4u1,n− u2,n = h2f1,n,
−u2,n−1− u1,n+ 4u2,n− u3,n = h2f2,n, ...
−un−1,n−1− un−2,n+ 4un−1,n− un,n = h2fn−1,n,
−un,n−1− un−1,n+ 4un,n = h2fn,n. Above equations can be represented as following matrix form:
−In A1
u:,n−1
u:,n
= h2f:,n.
師大
Therefore, (2) with boundary conditions is equivalent to a linear system Au = h2f with
A =
A1 −In
−In A1 . ..
. .. . .. −In
−In A1
∈ Rn2×n2, (6)
and
A1=
4 −1
−1 . .. ...
. .. ... −1
−1 4
, u =
u:,1 u:,2 ... u:,n
, f =
f:,1 f:,2 ... f:,n
.
師大
Question
How to compute the eigenpair (λ, x) of the matrix A which λ is closest to a target value?
Exercise
Prove that the vector f is an eigenvector of J = (4I − A)/4, also an eigenvector of A. Furthermore,
J f = µf and Af = 4(1 − µ)f with
µ = cos πh.
