Use of Chebyshev Polynomials to Construct L2 and Hilbert Space ...
Use of Chebyshev Polynomials to Construct L2 and Hilbert Space Inner Product to Solve Eigenvalue Problems With High Accuracy
A Project Submitted to the Faculty of the Department of Mathematics and Statistics University of Minnesota, Duluth by
Zhengfei Rui In partial fulfillment of the requirements for the degree of Master of Science in Applied and Computational Mathematics
Advisor: Steven A. Trogdon
June 24, 2014
Acknowledgments First and foremost, I would like to show my deepest gratitude to my supervisor, Dr. Steven A. Trogdon, a respectable, responsible and resourceful professor, who has provided me with valuable guidance in every stage of the writing of this project. Without his enlightening instruction, impressive kindness and patience, I could not have completed it. His keen and vigorous academic observation enlightens me not only in this paper work but also in my future study. Secondly, I shall extend my thanks to Dr. Zhuangyi Liu, Dr. Jiann-Shiou Yang and Dr. Debao Zhou for all their kindness and help. I would also like to thank all my teachers who have helped me to develop the fundamental and essential academic competence. Last but not least, I’ d like to thank my elder brother Zhenhua Rui for his encouragement and support. Zhengfei Rui
1
Abstract Chebyshev polynomials may be used to construct basis functions which satisfy fairly general boundary conditions. For the purpose of this report, these boundary conditions are associated with finding solutions to a system of coupled partial differential equations in space and time. These basis functions along with L2 and Hilbert space inner products may be used to spatially discretize the system of partial differential equations yielding a system of ordinary differential equations in time. By exploiting properties of Chebyshev polynomials the system of ordinary differential equations may be converted to an algebraic eigenvalue problem which may be solved using standard software. Using basis functions that satisfy boundary conditions is shown to give more accurate results for eigenvalues than when the Tau method, where basis function only approximately satisfy boundary conditions, is used.
2
Contents Acknowledgments
1
Abstract
2
List of Tables
4
List of Figures
5
List of Matlab Code
6
1 Introduction
7
2 Notation and Formulae
8
3 Beam Equations 3.1 Construction of a Chebyshev Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Spatially Discretized Equations Using a L2 Inner Product . . . . . . . . . . . . . . . 3.3 Spatially Discretized Equations Using a Hilbert Space Inner Product . . . . . . . . .
10 10 11 14
4 Numerical results 4.1 Decoupled problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Coupled problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21 21 21 25
A Appendix
27
B Matlab Code B.1 main . . . . . . . . . . . B.2 c . . . . . . . . . . . . . B.3 ChebyshevExpansion . . B.4 ChebyshevMultiplication B.5 GetClosestPoint . . . . . B.6 DrawSave1 . . . . . . . B.7 DrawSave2 . . . . . . . B.8 Tau1Beam . . . . . . . . B.9 Tau2Beam . . . . . . . . B.10 L2NB1Beam . . . . . . B.11 L2NB2Beam . . . . . . B.12 HSNB1Beam . . . . . . B.13 HSNB2Beam . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
3
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
29 29 31 31 31 32 32 33 34 35 36 37 38 40
List of Tables 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9
Decoupled result of Tau method . . . . . . . . . . . . . . . . . . . . . . . . . . . . Decoupled result of L2 inner product method based on Chebyshev basis . . . . . . Decoupled result of Hilbert Space inner product method based on Chebyshev basis Coupled result of Tau method, δ = 0 . . . . . . . . . . . . . . . . . . . . . . . . . Coupled result of Tau method, δ = 1 . . . . . . . . . . . . . . . . . . . . . . . . . Coupled result of L2 inner product method based on Chebyshev basis, δ = 0 . . . Coupled result of L2 inner product method based on Chebyshev basis, δ = 1 . . . Coupled result of Hilbert Space inner product method based on Chebyshev basis, δ Coupled result of Hilbert Space inner product method based on Chebyshev basis, δ
4
. .
22 22 22 . 23 . 23 . 24 . 24 = 0 24 = 1 25
List of Figures 4.1 4.2 4.3 4.4
Last Last Last Last
non-zero non-zero non-zero non-zero
eigenvalue, eigenvalue, eigenvalue, eigenvalue,
coupled coupled coupled coupled
problem problem problem problem
5
-
Tau method, δ = 0 . . . . . . . . . . Tau method, δ = 1 . . . . . . . . . L2 inner product method, δ = 1 . . Hilbert Space inner product method,
. . . δ
. 25 . 26 . 26 = 0 27
List of Matlab Code B.1 B.2 B.3 B.4 B.5 B.6 B.7 B.8 B.9 B.10 B.11 B.12 B.13
main.m . . . . . . . . . . . c.m . . . . . . . . . . . . . . ChebyshevExpansion.m . . ChebyshevMultiplication.m GetClosestPoint.m . . . . . DrawSave1.m . . . . . . . . DrawSave2.m . . . . . . . . Tau1Beam.m . . . . . . . . Tau2Beam.m . . . . . . . . L2NB1Beam.m . . . . . . . L2NB2Beam.m . . . . . . . HSNB1Beam.m . . . . . . . HSNB2Beam.m . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
6
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
29 31 31 31 32 32 33 34 35 36 37 38 40
1
Introduction
The spectral method is widely used to solve partial differential equations (PDEs). The method consists in approximating the solution to a PDE as a linear combination of a set of orthogonal basis functions. This approximation is called a spectral approximation. Once the approximation is substituted into the PDE, a projection procedure is employed to convert the original PDE into a system of equations for coefficients that multiply the basis functions in the spectral approximation. Chebyshev polynomials are often chosen as the basis functions used in implementing the spectral method. These polynomials are defined as Tn (x) = cos(n arccos x), −1 ≤ x ≤ 1, n = 0, 1, . . . and are orthogonal with respect to the inner product Z 1 uv √ dx hu, vi = 1 − x2 −1
(1.1)
(1.2)
In general, for a PDE of order 2k Lu = g, −1 ≤ x ≤ 1
(1.3)
with boundary conditions Bl,i u(−1) = 0, Br,i u(1) = 0, i = 1, . . . k
(1.4)
where L, Bl,i , Br,i are linear differential operators and g is a known function of x, we look for an approximate solution in the N + 1 dimensional Chebyshev function space VN = span{Tn (x), n = 0, . . . N } as u ∼ uN =
N X
an Tn (x)
(1.5)
(1.6)
n=0
and require that LuN − g be orthogonal to VN with respect to our inner product. This orthogonality or projection requirement means that hLuN − g, Ti (x)i = 0, i = 0, . . . N
(1.7)
which is equivalent to N X
n=0
an hLTn (x), Ti (x)i = hg, Ti (x)i, i = 0, . . . N
(1.8)
To satisfy the boundary conditions (1.4), we require that Bl,i uN (−1) = 0, Br,i uN (1) = 0, i = 1, . . . k
(1.9)
which is equivalent to N X
n=0
an Bl,i Tn (−1) = 0,
N X
an Br,i Tn (1) = 0, i = 1, . . . k
(1.10)
n=0
In equations (1.8) and (1.10), there are N + 1 + 2k equations with N + 1 unknowns. To make the number of equations and unknowns consistent, we apply the inner product in (1.8) for i = 0, . . . , N − 2k; the remaining 2k equations in (1.8) are dropped from the analysis.
7
2
Notation and Formulae For convenience, we introduce some notation. Suppose
f (x) =
N X
fi Ti (x) = f T TN
(2.1)
i=0
where
T0 (x) f0 f = ... , TN = ...
(2.2)
TN (x)
fN
and f T denotes the transpose of f . From now on, f (x) = f T TN is denoted as
V
N f f (x) −→
(2.3)
where VN is the Chebyshev function space defined in the Introduction (1.5). This notion means that f (x) may be expanded in the function space VN as the vector f . Furthermore, using properties of Chebyshev polynomials, it can be shown that V
N f⇒ f (x) −→
dk VN M kf , k ≥ 1 f (x) −→ k dx
(2.4)
where the matrix M is of order (N + 1) × (N + 1) and is defined in the Appendix (A.5). Now let g0 VN .. g(x) −→ g = . .
(2.5)
gN
By using the multiplication law,
Ti (x)Tj (x) =
1 1 Ti+j (x) + T|i−j| (x) 2 2
(2.6)
for the Chebyshev polynomials, the product f (x)g(x) may be approximated as f (x)g(x)
= = = =
N X
fi Ti (x)
i=0 N N X X
i=0 j=0 N N X X
N X
fi gj Ti (x)Tj (x) fi g j
i=0 j=0 2N X
1 2
gj Tj (x)
j=0
1 1 Ti+j (x) + T|i−j| (x) 2 2
S1,i Ti (x) +
i=0
N
1X S2,i Ti (x) 2 i=0
where S1,i =
X
fp g p ,
p+q =i 0 ≤ p, q ≤ N
8
0 ≤ i ≤ 2N
(2.7)
and X
S2,i =
0≤i≤N
fp g p ,
|p − q| = i 0 6 p, q 6 N
(2.8)
From (2.7) we have that S1,0 S1,1 S1,2 .. .
= = =
f0 g 0 f0 g 1 + f1 g 0 f0 g 2 + f1 g 1 + f2 g 0
S1,N
=
f0 gN + f1 gN −1 + · · · + fN g0
and thus
S1,0 f0 S1,1 f1 .. = .. . . S1,N
f0
fN
fN −1
···
Likewise from (2.8) we have that S2,0 S2,1 .. . S2,N
= f0 g 0 + f1 g 1 + · · · + fN g N = f0 g1 + f1 g2 + · · · + fN −1 gN + f1 g0 + f2 g1 + · · · + fN gN −1
(2.9)
(2.10) (2.11) (2.12) (2.13)
= f0 g N + fN g 0
and thus that S2,0 f0 S2,1 f1 .. = .. . . S2,N
g 0 .. . ≡ W1 g. .. . f0 gN
fN
f1 f0 + f2
f2 f1 + f3
··· ···
fN −1 fN −2 + fN
0
···
···
0
Hence we get that
f (x)g(x)
= = = =
g 0 .. fN −1 . ≡ W2 g. .. . f0 gN fN
(2.14)
2N N 1X 1X S1,i Ti (x) + S2,i Ti (x) 2 i=0 2 i=0 N N 2N 1X 1 X 1X S1,i Ti (x) + S2,i Ti (x) + S1,i Ti (x) 2 i=0 2 i=0 2 i=N +1 1 1 1 T T T [W1 g] TN + [W2 g] TN + [W3 g] TN∗ 2 2 2 1 1 T T [(W1 + W2 )g] TN + [W3 g] TN∗ 2 2
where TN∗
TN +1 (x) .. = . .
(2.15)
T2N (x)
If we let Wf = 21 (W1 + W2 ), and Wf∗ = 12 W3 , then we get that T T f (x)g(x) = [Wf g] TN + Wf∗ g TN∗
(2.16)
where Wf∗ is not given explicitly because it consists of only terms of order higher than N , which become 0 after taking the inner product with Tk (x), k = 0, 1, . . . N. 9
3
Beam Equations
We are given the following system of partial differential equations describing the motion of two coupled beams ∂ 4 u1 ∂u1 ∂ 2 u1 = −p1 − k(x) (u1 − u2 ) − δk(x) , 2 ∂t ∂x4 ∂t ∂ 2 u2 ∂ 4 u2 ∂u2 = −p + k(x) (u1 − u2 ) − δk(x) . 2 ∂t2 ∂x4 ∂t
(3.1a) (3.1b)
with boundary conditions
2
u1 (t, −1) =u1 (t, 1) = 0,
(3.2a)
2
∂ u1 ∂ u1 (t, −1) = (t, 1) = 0, 2 ∂x ∂x2 u2 (t, −1) =u2 (t, 1) = 0,
(3.2b) (3.2c)
2
2
∂ u2 ∂ u2 (t, −1) = (t, 1) = 0. ∂x2 ∂x2
(3.2d)
where p1 > 0, p2 > 0, and δ ≥ 0 are given parameters and k(x) is a known function. ∂u2 1 If we introduce the substitutions, v1 = ∂u ∂t and v2 = ∂t , then the original equations (3.1a) and (3.1b) become ∂u1 ∂t ∂v1 ∂t ∂u2 ∂t ∂v2 ∂t
(3.3a)
= v1 , 4
= −p1
∂ u1 − k(x) (u1 − u2 ) − δk(x)v1 , ∂x4
(3.3c)
= v2 , = −p2
(3.3b)
∂ 4 u2 + k(x) (u1 − u2 ) − δk(x)v2 . ∂x4
(3.3d)
Next, in section 3.1 we will construct basis functions that satisfy the boundary conditions using Chebyshev polynomials. Then, in sections 3.2 and 3.3 we use two different projection methods to set up eigenvalue (linear algebra) problems which can be solved, and in section 4 we compare the numerical results these two different projection methods give us.
3.1
Construction of a Chebyshev Basis
In this section we construct a set of functions that span a vector space GN , where each basis element of GN satisfies boundary conditions (3.2a)-(3.2d). We then expand u1 ,v1 ,u2 ,v2 in GN as they should satisfy the associated boundary conditions, and solve (3.3a)-(3.3d) with boundary conditions (3.2a)-(3.2d). Let GN = span{Gn (x), n = 0, ...N } where Gn = Tn (x) + a1 Tn+1 (x) + a2 Tn+2 (x) + a3 Tn+3 (x) + a4 Tn+4 (x) = Tn (x) +
4 X
ai Tn+i (x)
i=1
10
Plugging the above Gn (x) into (3.2a) and (3.2b), and using formula (A.2) from Appendix A, we get
Gn (−1) = 0 ⇒ 1 + Gn (1) = 0 ⇒ 1 +
4 X
(−1)i ai = 0,
(3.5a)
ai = 0,
(3.5b)
i=1
4 X i=1
4
2 2 d2 n2 (n2 − 1) X i (n + i) [(n + i) − 1]ai G (−1) = 0 ⇒ (−1) + = 0, n dx2 3 3 i=1
(3.5c)
4
n2 (n2 − 1) X (n + i)2 [(n + i)2 − 1]ai d2 (−1)i Gn (−1) = 0 ⇒ + = 0. 2 dx 3 3 i=1
(3.5d)
These are a set of algebraic equations with four unknowns, ai , i = 1, ...4. We use Mathematica to solve for these unknowns and get that a1 (n) = a3 (n) = 0, a2 (n) = −
2(n + 2)(2n2 + 8n + 15) (n + 1)(2n2 + 4n + 3) , a (n) = . 4 (n + 3)(2n2 + 12n + 19) (n + 3)(2n2 + 12n + 19)
(3.6)
Therefore, Gn (x) = Tn (x) + a2 TN +2 (x) + a4 Tn+4 (x).
(3.7)
It can be shown that G0 (x), G1 (x), ...GN (x) are linear independent, and thus span the subspace GN of VN +4 . It is easy to see that, every g(x) ∈ GN is representable in VN +4 . In fact, we have that VN +4
G
N g ⇒ g(x) −→ Gg g(x) −→
(3.8)
where
1
0 a2 (0) 0 a4 (0) G= 0 . .. . .. . .. 0
0
0
1
0
0
1
a2 (1) 0 a4 (1)
0 a2 (2) 0
.. . .. . 0
a4 (2) 0 .. . 0
···
0 .. . .. . .. .
0 . 1 0 a2 (N ) 0
(3.9)
a4 (N )
Thus G has order (N + 5) × (N + 1) and represents the linear transformation from GN to VN +4 .
3.2
Spatially Discretized Equations Using a L2 Inner Product
Let’s map the u1 , u2 , v1 , v2 to the VN +4 , so that we can apply the product formula (2.16) and inner product relation (A.3) to them with Tn (x), n = 0, . . . , N . We also expand k(x) in VN +4 since it will
11
multiply functions from VN +4 instead of GN . The expansions are as follows: V N +4 a0 (t) u1 −→ Ga(t) GN VN +4 1 a(t) = ... ⇒ ∂u u1 −→ ˙ ∂t −→ Ga(t) ∂ 4 u1 VN +4 4 aN (t) ∂x4 −→ M Ga(t), VN +4 b0 (t), u2 −→ Gb(t) GN VN +4 2 ˙ b(t) = ... ⇒ ∂u u2 −→ ∂t −→ Gb(t) 4 ∂ u2 VN +4 4 bN (t) ∂x4 −→ M Gb(t) c0 (t), N +4 v V−→ Gc(t) GN 1 c(t) = ... ⇒ v1 −→ N +4 ∂v1 V−→ Gc(t) ˙ ∂t cN (t) d0 (t) N +4 v V−→ Gd(t) GN .. 2 v2 −→ d(t) = . ⇒ N +4 ∂v2 V−→ ˙ Gd(t), ∂t dN (t) k0 . VN +4 k(x) −→ k = ... , kN +4
where the superposed dot denotes time differentiation. Upon substituting from the above expansions into (3.3a), we get that N X
a˙i (t)Gi (x) =
N X
ci (t)Gi (x)
i=0
i=0
and since the Gi (x), i = 0, . . . , N are linearly independent, we deduce that a˙ i (t) = ci (t), i = 0, . . . , N. This can be written in matrix-vector form as IN +1 a(t) ˙ = IN +1 c(t)
(3.11)
where IN +1 is the identity matrix of order N +1. This is the spatially discretized equation associated with (3.3a). Likewise, we can also get the spatially discretized equation associated with (3.3c) as ˙ = IN +1 d(t). IN +1 b(t)
(3.12)
The remaining equations to be discretized are equations (3.3b) and (3.3d). Using (2.3) for functions expanded in the function space, VN +4 we get that k(x)u1 (x) = [Wk Ga(t)]T TN +4 + [Wk∗ Ga(t)]T TN∗ +4 , k(x)u2 (x) = [Wk Gb(t)]T TN +4 + [Wk∗ Gb(t)]T TN∗ +4 , k(x)v1 (x) = [Wk Gc(t)]T TN +4 + [Wk∗ Gc(t)]T TN∗ +4 , k(x)v2 (x) = [Wk Gd(t)]T TN +4 + [Wk∗ Gd(t)]T TN∗ +4 . With the function space expansion of u1 , u2 and v1 in VN +4 , we then obtain from (3.3b) that [Gc(t)] ˙ T TN +4 = [−p1 M 4 Ga(t)]T TN +4 − [Wk Ga(t)]T TN +4 − [Wk∗ Ga(t)]T TN∗ +4 ,
+[Wk Gb(t)]T TN +4 + [Wk∗ Gb(t)]T TN∗ +4 − δ[Wk Gc(t)]T TN +4 − δ[Wk∗ Gc(t)]T TN∗ +4 . 12
After taking the inner product of these equations with Tj (x), j = 0, . . . N + 4 we are left with N + 5 equations in N + 1 unknowns. To make the system of equations consistent, we keep the first N + 1 equations and drop the last 4 equations. In fact, if k(x) = 0, i.e. the two beams are independent of each other, the right hand side of the last 4 equations is identically zero which makes the system of equation inconsistent unless those equations are dropped. Upon taking the inner product with Tj (x), j = 0, . . . N , we get ˙ = [−p1 M 4 G]j a(t) − [Wk G]j a(t) + [Wk G]j b(t) − δ[Wk G]j c(t), j = 0, . . . N. [G]j c(t) This can be written in matrix-vector form as ˙ = A1 a(t) − Ba(t) + Bb(t) − δBc(t). G1 c(t)
(3.13)
Equations (3.13) are the spatially discretized equations corresponding to (3.3b), where G1 , of order(N + 1) × (N + 1) is the submatrix obtained from the first N + 1 rows of G. And A1 , of order(N + 1) × (N + 1) is the submatrix obtained from the first N + 1 rows of −p1 M 4 G. Also B, of order(N + 1) × (N + 1) is the submatrix obtained from the first N + 1 rows of Wk G. In the same way, we get ˙ = [−p2 M 4 G]j b(t) + [Wk G]j a(t) − [Wk G]j b(t) − δ[Wk G]j d(t), j = 0, . . . N [G]j d(t) This can also be written in matrix-vector form as, ˙ = A2 b(t) + Ba(t) − Bb(t) − δBd(t). G1 d(t)
(3.14)
where A2 , of order(N +1)×(N +1) is the submatrix obtained from the first N +1 rows of −p2 M 4 G. If we let a(t) c(t) (3.15) z(t) = b(t) d(t)
then (3.11)-(3.14) may be combined as the following
IN +1 0 0 0
0 0 G1 0
0 IN +1 0 0
0 0 0 0 z(t) ˙ = A1 − B 0 B G1
IN +1 0 −δB 0
0 0 B A2 − B
which may more compactly be written as
˙ Qz(t) = Dz(t).
0
IN +1 z(t) 0 −δB (3.16)
The above equations are the spatially discretized system of equations associated with using basis functions which satisfying the boundary conditions. The eigenvalue problem for these equations is ˆ where obtained by setting z(t) = eλt z, ˆ a ˆb0 a ˆ0 cˆ0 dˆ0 cˆ . ˆ = .. , cˆ = .. , dˆ = .. . ˆ= zˆ = (3.17) .. , b . . . ˆ , a b ˆ ˆ a ˆN cˆN bN dN dˆ Therefore we obtain
λQeλt zˆ = Deλt zˆ ⇒ λzˆ = Q−1 D zˆ which is the eigenvalue problem corresponding to the equations (3.16). 13
(3.18)
3.3
Spatially Discretized Equations Using a Hilbert Space Inner Product
The Hilbert space inner product that we use in this section requires that we calculate the L2 inner product, hD 2 Gi , D 2 Gj iL2 . From equations (2.4) and (A.5) it is fairly straightforward to show that D 2 Gi becomes quite large as the index, i becomes large and thus elements in the matrix representing the above inner product can become quite large. It is the u1 and u2 in equations (3.3a)(3.3d) that are needed in calculating the elements of hD 2 Gi , D 2 Gj iL2 . Thus, in order to insure that all variables have the same relative order of magnitude, we introduce a scaling parameter, ri to normalize the coefficients of the vectors and let ai (t) ri u1 v 1 GN ci (t) (3.19) z= , u2 −→ bi (t) v2 ri di (t) where the following series expansions
u1 =
N X
ai (t)
i=0
v1 =
N X
Gi (x) , ri
ci (t)Gi (x),
i=0
u2 =
N X
bi (t)
i=0
v2 =
N X
Gi (x) , ri
di (t)Gi (x)
i=0
are used and where ri , the scaling parameter is defined in equation (A.6) in Appendix A. Upon substituting the above expressions for u1 , v1 , u2 , v2 into (3.3a)-(3.3d), we get that N X
N
Gi (x) X a˙i (t) ci (t)Gi (x), = ri i=0 i=0 N X i=0
N X
c˙i (t)Gi (x) = −p1
Gi (x) = b˙i (t) ri i=0 N X i=0
N X
N X
ai (t)
i=0
(3.20)
N N X X Gi (x) D 4 Gi (x) (ai (t) − bi (t)) ci (t)Gi (x), − k(x) − δk(x) ri ri i=0 i=0
(3.21)
(3.22)
di (t)Gi (x),
i=0
d˙i (t)Gi (x) = −p2
N X i=0
bi (t)
N N X X D 4 Gi (x) Gi (x) (ai (t) − bi (t)) di (t)Gi (x). + k(x) − δk(x) ri ri i=0 i=0
(3.23)
14
Now let the left hand side (LHS) of the above equations be N X
Gi (x) a˙ i (t) ri i=0 N X c˙i (t)Gi (x) i=0 . zl = N X G (x) i b˙ i (t) r i i=0 X N d˙i (t)Gi (x)
(3.24)
i=0
from which we may define the following four vectors as a basis for representing the elements of zl
zˆ(1)
Gj (x) 0 0 0 rj 0 (2) Gj (x) (3) (4) 0 = 0 , zˆ = 0 , zˆ = Gj (x) , zˆ = 0 . 0 rj 0 Gj (x) 0 0
(3.25)
In order to obtain the best approximation, we will utilize a suitably defined Hilbert space inner product [1, 4] to spatially discretize the system of equations given above. We provide the following theorems that relate to the system of equations in section 3 as justification for the use of a Hilbert space inner product: Theorem 3.1 The operator A generates a semigroup, eAt , on the Hilbert space H. Therefore, the system is well-posed, i.e, for any Z0 ∈ H, the unique solution to the system is Z(t) = eAt Z0 , t > 0. [4] Theorem 3.2 The semigroup eAt is exponentially stable. [4] where H = U1 × V 1 × U2 × V 2 , U1 = H01 (−1, 1) ∩ H 2 (−1, 1) , V1 = L2 (−1, 1),
and
U2 = H01 (−1, 1) ∩ H 2 (−1, 1) , V2 = L2 (−1, 1), 1
f 2 (x) √ L (−1, 1) = f (x) | hf, f iL2 ≡ dx < ∞ , 1 − x2 −1 n o H k (−1, 1) = f (x) | f (x), . . . , f (k−1) (x) ∈ L2 (−1, 1) , H01 (−1, 1) = f (x) ∈ H 1 (−1, 1) | f (−1) = f (1) = 0 . 2
Z
Thus, relative the beam equations in section 3, for any z = (u1 , v1 , u2 , v2 )T ∈ H and any ˆ z = (ˆ u1 , vˆ1 , u ˆ2 , vˆ2 )T ∈ H, we may define the following inner product:
ˆ H = (u1 , v1 , u2 , v2 )T , (ˆ hz, zi u1 , vˆ1 , u ˆ2 , vˆ2 )T H ≡
p1 hD 2 u1 , D 2 u ˆ1 iL2 + p2 hD 2 u2 , D 2 u ˆ2 iL2 + hv1 , vˆ1 iL2 + hv2 , vˆ2 iL2 +hk(x)(u1 − u2 ), (uˆ1 − u ˆ2 )iL2 .
15
We now compute hzl , zˆ(1) iH , hzl , zˆ(2) iH , hzl , zˆ(3) iH and hzl , zˆ(4) iH , where the simplified notation, h·, ·i ≡ h·, ·iL2 is employed: N N a˙ i − b˙i X X a ˙ i hD 2 Gi , D 2 Gj i + hk(x)Gi , Gj i, (3.26) hzl , zˆ(1) iH = p1 r r ri rj i j i=0 i=0 hzl , zˆ(2) iH = hzl , zˆ(3) iH
N X i=0
(3.27)
c˙i hGi , Gj i,
N N a˙ i − b˙i X X ˙bi = p2 hD 2 Gi , D 2 Gj i − hk(x)Gi , Gj i, r r ri rj i j i=0 i=0
hzl , zˆ(4) iH =
N X i=0
d˙i hGi , Gj i.
(3.28)
(3.29)
Now let
N X
ci Gi (x) i=0 N N N 4 X X X D Gi (x) Gi (x) ai (ai − bi ) ci Gi (x) − k(x) − δk(x) −p1 r r i i i=0 i=0 i=0 zr = N X di Gi (x) i=0 N N N X X X D 4 Gi (x) Gi (x) −p2 (ai − bi ) di Gi (x) bi + k(x) − δk(x) ri ri i=0 i=0 i=0
(3.30)
be the right hand side (RHS) of our system of equation and compute hzr , zˆ(1) iH , hzr , zˆ(2) iH , hzr , zˆ(3) iH and hzr , zˆ(4) iH , N N c˙i − d˙i X X c ˙ i hzr , zˆ(1) iH = p1 hD 2 Gi , D 2 Gj i + hGi , Gj i, (3.31) r rj j i=0 i=0 ˙ N N N a ˙ − b X X X i i a˙ i hzr , zˆ(2) iH = −p1 ci hk(x)Gi , Gj i, (3.32) hD 4 Gi , Gj i − hk(x)Gi , Gj i − δ r ri i=0 i i=0 i=0 ˙i N ˙ N c ˙ − d X X i di hzr , zˆ(3) iH = p2 hD 2 Gi , D 2 Gj i − hG, Gj i, (3.33) r rj j i=0 i=0 N N N ˙ a˙ i − b˙i X X X b i di hk(x)Gi , Gj i. (3.34) hD 4 Gi , Gj i + hk(x)Gi , Gj i − δ hzr , zˆ(4) iH = −p2 r ri i=0 i=0 i=0 i
16
The inner product operation applied to the LHS should be equal to that applied to the RHS, thus ˙ N N a˙ i (t) − bi (t) X X a˙ i (t) 2 p1 hD Gi , D 2 Gj i + hk(x)Gi , Gj i r r ri rj i j i=0 i=0 ˙i (t) N N c ˙ (t) − d X X i c˙i (t) 2 hD Gi , D 2 Gj i + hGi , Gj i, (3.35a) = p1 r rj j i=0 i=0 ˙i (t) N N N a ˙ (t) − b X X X i a˙ i (t) 4 c˙i (t)hGi , Gj i = −p1 hD Gi , Gj i − hk(x)Gi , Gj i ri ri i=0 i=0 i=0 N X
−δ p2
N ˙ X bi (t) i=0
ri rj
hD 2 Gi , D 2 Gj i − = p2
i=0
d˙i (t)hGi , Gj i = −p2
ri rj
i=0
N ˙ X di (t)
rj
i=0
N X
N a˙ i (t) − b˙i (t) X
ri
hD 4 Gi , Gj i +
(3.35b)
ci (t)hk(x)Gi , Gj i,
hk(x)Gi , Gj i
hD 2 Gi , D 2 Gj i −
N ˙ X bi (t) i=0
i=0
N c˙i (t) − d˙i (t) X rj
i=0
N a˙ i (t) − b˙i (t) X ri
i=0
−δ
N X i=0
(3.35c)
hGi , Gj i,
hk(x)Gi , Gj i (3.35d)
di (t)hk(x)Gi , Gj i.
Here we will perform some calculations to simplify items in the above equations. For every g(x) ∈ GN , g(x) can expressed in terms of either Gi or Ti as follows: g(x) =
N X
gi Gi (x) =
N +4 X
(3.36)
gˆi Ti (x), i = 0, . . . , N + 4.
i=0
i=0
From (3.8), we know that gˆ = Gg so that gˆi =
N X
(3.37)
Gi,k gk , i = 0. . . . , N + 4,
k=0
and g(x) =
N N +4 X X
Gi,k gk Ti (x) =
i=0 k=0
N X
k=0
N +4 X
!
(Gi,k Ti (x)) gk =
i=0
N X i=0
N +4 X
!
(Gk,i Tk (x)) gi .
k=0
(3.38)
Comparing the above to (3.36), it is easy to see that Gi (x) =
N +4 X
(3.39)
Gk,i Tk (x), i = 0, . . . , N.
k=0
So if we define GN = GT TN +4 , where G is (N + 5) × (N + 1) and TN +4 is (N + 5) × 1, then from (3.36) g(x) = g T GN . Likewise, for the pth derivative g
(p)
(x) =
N X i=0
(p) gi Gi (x)
=
N +4 X i=0
17
(p) gˆi Ti (x)
=
N +4 X i=0
(p)
gˆi Ti (x)
(3.40)
where from equation (2.4) (p)
gˆi
=
N +4 X
(M p )i,k gˆk .
(3.41)
k=0
So we can express g (p) (x) as follows: +4 N +4 N X X
g (p) (x) =
(M p )i,k gˆk Ti (x)
i=0 k=0
+4 N +4 N X X
=
(M p )i,k
=
+4 N +4 N X X
m=0 N X
=
Gk,m gm Ti (x)
m=0
i=0 k=0
N X
N X
!
(M p )i,k Gk,m Ti (x) gm
i=0 k=0
G(p) m (x)gm .
m=0
Now if we switch the indices m and i, we get that g (p) (x)
=
N P
i=0 (p)
Gi (x) =
+4 NP +4 NP m=0 k=0
(p)
Gi (x)gi , where
(p)
(M p )m,k Gk,m Tm (x), i = 0, ..., N . Therefore if we define GN = (M p G)T TN +4
we may write (3.40) as (p)
g (p) (x) = g T GN .
(3.42)
Now we simplify the other inner products that appear in equations (3.35a) - (3.35d). *N +4 + N +4 N +4 N +4 X X X X Gk,i Gm,j hTk (x), Tm (x)i hGi , Gj i = Gm,j Tm (x) = Gk,i Tk (x), k=0 m=0
m=0
k=0
Since hTk (x), Tm (x)i = αk δk,m , then when k = m, hTk (x), Tm (x)i = αk . So that hGi , Gj i = ˆ k,i = √αk Gk,i , then If we define G hGi , Gj i =
N +4 X
N +4 X
(3.43)
αk Gk,i Gk,j .
k=0
ˆ k,i G ˆ k,j = G
k=0
N +4 X
ˆ T )i,k G ˆ k,j . (G
(3.44)
k=0
In the same way, 2
2
hD Gi , D Gj i = =
*N +4 X
(M
2
G)T i,k Tk (x),
(M
2
G)T j,m Tm (x)
m=0
k=0
+4 N +4 N X X
N +4 X
+
(3.45)
2 T (M 2 G)T i,k (M G)j,m hTk (x), Tm (x)i
(3.46)
2 T (M 2 G)T i,k (M G)j,m αk δk,m .
(3.47)
k=0 m=0
=
+4 N +4 N X X
k=0 m=0
18
ˆ 2 = √αk M 2 , then and if we define M k,i k,i 2
N +4 X
2
hD Gi , D Gj i =
ˆ 2 G)T (M ˆ 2 G)T = (M i,k j,k
N +4 X
ˆ 2 G)T (M ˆ 2 G)k,j . (M i,k
(3.48)
k=0
k=0
ˆ 4 = √αk M 4 , we get that Now defining M k,i k,i *N +4 + N +4 N +4 X X X ˆ 4 G)T G ˆ hD 4 Gi , Gj i = (M (G)T = (M 4 G)T i,k k,j . j,m Tm (x) i,k Tk (x), m=0
k=0
(3.49)
k=0
Lastly, the inner product hk(x)Gi , Gj i will be calculated, k(x)Gi (x) = k(x)GT TN +4 =
N +4 X
k(x)Gk,i Tk (x).
(3.50)
k=0
If we substitute k(x) =
PN +4
m=0
km Tm (x) into the above equation then we get that
k(x)Gi (x) =
+4 N +4 N X X
km Gk,i Tm (x)Tk (x).
(3.51)
k=0 m=0
Since Tm (x)Tk (x) = 21 Tm+k (x) + 12 T|m−k| (x), then we may write the above equation as T
T
k(x)Gi (x) = [Wk gi ] TN +4 + [Wk∗ gi ] TN∗ +4 =
N +4 X
[Wk gi ]m Tm (x) +
N +4 X
∗ [Wk∗ gi ]m Tm+N +5 (x)
(3.52)
m=0
m=0
PN +4 where [Wk gi ]m = n=0 (Wk )m,n (gi )n . We will neglect the Wk∗ quantities since these terms will contribute nothing when taking the inner product with Gj (x). Since (gi )n = Gn,i , then k(x)Gi (x) =
+4 N +4 N X X
(Wk )m,n Gn,i Tm (x)
(3.53)
m=0 n=0
Upon taking the inner product of the above equation with Gj (x) we get that + *N +4 N +4 N +4 X X X Gp,j Tp (x) hk(x)Gi , Gj i = (Wk )m,n Gn,i Tm (x), p=0
m=0 n=0
=
+4 +4 N N +4 N X X X m=0 n=0 p=0
=
+4 N +4 N X X
(Wk )m,n Gn,i Gp,j hTm (x), Tp (x)i
αm (Wk )m,n Gn,i Gm,j .
(3.54)
m=0 n=0
ˆ k )m,n = √αm (Wk )m,n with G ˆ m,j = √αm Gm,j , the above equation may And thus after defining (W be expressed as hk(x)Gi , Gj i =
N +4 X
ˆ k G)T G ˆ (W i,m m,j .
(3.55)
m=0
Thus the left-hand side of equations (3.35a)-(3.35d) may be written in block matrix form as ˙ a(t) M11 0 M13 0 ˙ 0 M22 0 0 c(t) (3.56) M z˙ = ˙ , b(t) M31 0 M33 0 ˙ 0 0 0 M44 d(t) 19
where the matrix blocks are given by hD 2 Gi , D 2 Gj i hk(x)Gi , Gj i + , ri rj ri rj hk(x)Gi , Gj i , = ri rj
(M11 )j,i = p1 (M13 )j,i
(M22 )j,i = hGi , Gj i,
hk(x)Gi , Gj i = (M13 )j,i ri rj hD 2 Gi , D 2 Gj i hk(x)Gi , Gj i + , = p2 ri rj ri rj
(M31 )j,i = (M33 )j,i
(M44 )j,i = hGi , Gj i = (M22 )j,i . and the right-hand side of equations (3.35a)-(3.35d) 0 A12 0 A21 A22 A23 Az = 0 A32 0 A41 0 A43 where the matrix blocks are given by
may be written in block matrix form as a(t) A14 0 c(t) , (3.57) b(t) A34 d(t) A44
hD 2 Gi , D 2 Gj i hGi , Gj i + , rj rj hGi , Gj i =− , rj hD 4 Gi , Gj i hk(x)Gi , Gj i − , = −p1 ri ri = −δhk(x)Gi , Gj i,
(A12 )j,i = p1 (A14 )j,i (A21 )j,i (A22 )j,i
hk(x)Gi , Gj i , ri hGi , Gj i = (A14 )j,i , =− rj
(A23 )j,i = (A32 )j,i
hD 2 Gi , D 2 Gj i hGi , Gj i + , rj rj hk(x)Gi , Gj i = = (A23 )j,i , ri hD 4 Gi , Gj i hk(x)Gi , Gj i − , = −p2 ri ri = −δhk(x)Gi , Gj i = (A22 )j,i .
(A34 )j,i = p2 (A41 )j,i (A43 )j,i (A44 )j,i
Thus upon equating equations (3.56) and (3.57) we obtain the ential equations ˙ a(t) 0 M11 0 M13 0 0 ˙ c(t) A M 0 0 22 21 M z˙ = Az ⇒ ˙ = 0 M31 0 M33 0 b(t) ˙ A41 0 0 0 M44 d(t) with corresponding algebraic eigenvalue problem
ˆ λzˆ = M −1 Az. 20
following system of ordinary differA12 A22 A32 0
0 A23 0 A43
a(t) A14 0 c(t) A34 b(t) A44 d(t)
(3.58)
(3.59)
4
Numerical results
All the numerical experiments for equations (3.18) and (3.59) were performed with Matlab 2012b(64bit) on a Windows 7 system. According to the Matlab documentation, the command eig(D, Q) provides higher accuracy than eig(Q−1 D), and eig(Q−1 D) provides higher accuracy than does eig(Q\D).
4.1
Decoupled problem
When k(x) = 0, then the two beams equations (3.1a)-(3.1b) will be simplified to be two beam equations with exactly the same boundary conditions. In this part, we will consider the simplest case where p1 = p2 = 1, k(x) = 0, then (3.3a)-(3.3b) become ∂u1 = v1 ∂t ∂v1 ∂ 4 u1 =− ∂t ∂x4
(4.1a) (4.1b)
with boundary conditions (3.2a)-(3.2b), and ∂u2 = v2 ∂t ∂v2 ∂ 4 u2 =− ∂t ∂x4
(4.2a) (4.2b)
with boundary conditions (3.2c)-(3.2d). The system for u1 and u2 are thus decoupled into two beam equations with identical boundary conditions. In the following results from simulations, we calculate eigenvalues for N = 20 to N = 200 and only focus on those eigenvalues with magnitude less than 300. We present, in tabular form, part of the simulation results where N is a multiple of 20, starting with N = 20 up to N = 200. For a given N , we define that eigenvalue having the smallest magnitude as the first non-zero eigenvalue and that eigenvalue having magnitude closest to 300 as the last non-zero eigenvalue. Tables (4.1)-(4.3) are for the decoupled problem and we list the first and last non-zero eigenvalue as given by the Tau method [2] and when using L2 and Hilbert Space inner products to spatially discretize our system of partial differential equations. When using the Tau method, the variables u1 , u2 , v1 , v2 are expanded in the VN space [3] while for the L2 and Hilbert Space projection methods the variables are expanded in the VN +4 space. We note that with the L2 and Hilbert Space projection methods, the basis functions satisfy the boundary condition exactly while with the Tau method boundary conditions are only satisfied approximately. From Tables (4.1)-(4.3) we see that, regardless of whether we focus on the first non-zero eigenvalue or the last non-zero eigenvalue, the number of significant digits in the imaginary part of the eigenvalues given by the Tau method is only 5. But from results where the L2 and Hilbert Space inner products were used to spatially discretize our system of partial differential equation, we have at least 8 significant digits. The real part of the eigenvalues should be identically zero. Numerical results when using the L2 inner product as a projection method returned this exact value while with the Hilbert Space inner product method real parts were on the order of zero, to machine accuracy. The Tau method returned the poorest results for the real part of the eigenvalues.
4.2
Coupled problem 1 2
1 2
In this section, we choose k(x) = 0.1(e−16(x− 2 ) + e−16(x+ 2 ) ), p1 = p2 = 1, with the damping parameter δ = 0 or δ = 1. Tables (4.4)-(4.9) are for the coupled problem and we list the first and last non-zero eigenvalue as given by the Tau method and when using L2 and Hilbert Space inner
21
Table 4.1: Decoupled result of Tau method first non-zero eigenvalue last non-zero eigenvalue -6.705598014098E-14+2.467401094780 i 8.847124421951E-13+2.986688399691E+02 i 1.568884491585E-11+2.467402097159 i 1.761809974169E-11+2.985555331487E+02 i 7.165028143492E-11+2.467423720752 i -2.349985916899E-10+2.985555322398E+02 i -1.503783152077E-10+2.467613460975 i 7.668656776916E-10+2.985555313579E+02 i -2.032218010365E-10+2.466537807889 i -8.733368354941E-10+2.985555333780E+02 i 4.437809932871E-10+2.462436402624 i -3.840691884082E-09+2.985555342675E+02 i 3.377024176437E-09+2.476829981371 i 5.128603730863E-09+2.985556439165E+02 i 8.018927187781E-10+2.514746929390 i -7.461088354724E-09+2.985539417344E+02 i -1.022905982853E-09+2.131379055191 i -1.732687196718E-08+2.985547721941E+02 i -6.756858552262E-09+2.337543707013 i -2.889222089258E-08+2.985567659389E+02 i
N 20 40 60 80 100 120 140 160 180 200
Table 4.2: Decoupled result of L2 inner product method based on Chebyshev basis N first non-zero eigenvalue last non-zero eigenvalue 20 0+2.467401100269 i 0+2.986688399692E+02 i 40 0+2.467401100246 i 0+2.985555331330E+02 i 60 0+2.467401100172 i 0+2.985555331330E+02 i 80 0+2.467401100228 i 0+2.985555331330E+02 i 100 0+2.467401100475 i 0+2.985555331330E+02 i 120 0+2.467401100171 i 0+2.985555331330E+02 i 140 0+2.467401099878 i 0+2.985555331330E+02 i 160 0+2.467401099521 i 0+2.985555331330E+02 i 180 0+2.467401100908 i 0+2.985555320838E+02 i 200 0+2.467401100218 i 0+2.985555347055E+02 i
Table N 20 40 60 80 100 120 140 160 180 200
4.3: Decoupled result of Hilbert Space inner product method based on Chebyshev basis first non-zero eigenvalue last non-zero eigenvalue -1.306999380261E-14+2.467401100174 i -9.98326702631E-11+2.991678018448E+02 i -1.343450278756E-13+2.467401101420 i 1.169260017025E-11+2.985555331311E+02 i -7.554568084249E-14+2.467401072860 i 1.196110695876E-10+2.985555331703E+02 i 7.987145869450E-14+2.467401172821 i -1.103474043019E-10+2.985555330376E+02 i 2.083246001098E-15+2.467401191054 i 1.925940433517E-11+2.985555330213E+02 i -1.620447780703E-14+2.467401361637 i -8.172333694982E-11+2.985555328481E+02 i 1.418744016787E-12+2.467398981886 i -5.641281441202E-11+2.985555356969E+02 i 6.989363144672E-11+2.467403232593 i -1.007735885850E-10+2.985555305601E+02 i 9.026536089681E-10+2.467398107616 i -6.140341590391E-09+2.985555367307E+02 i 6.478229430914E-11+2.467397967167 i 7.649350863156E-10+2.985555373950E+02 i
22
products to spatially discretize our system of partial differential equations with δ = 0 and δ = 1. From Tables (4.4)-(4.9) we see that, regardless of whether we focus on the first non-zero eigenvalue or the last non-zero eigenvalue, the number of significant digits in the imaginary part of the eigenvalues given by the Tau method is only 5. But from results where the L2 and Hilbert Space inner products were used to spatially discretize our system of partial differential equation, we have at least 7 significant digits. The real part of the eigenvalues should be identically zero when δ = 0 and should be strictly negative when δ = 1. When δ = 0 the Hilbert Space projection method returns eigenvalues with real part identically zero and the L2 projection method returns eigenvalues with real part on the order of zero, to machine accuracy. The Tau method returns the poorest approximation to the real part of the eigenvalues and exhibits some stability problems relative to convergence of the imaginary part of the eigenvalues. When δ = 1 both the L2 and Hilbert Space projection method return eigenvalues with real part that have the desired property of being strictly negative and it would appear that the L2 projection method returns eigenvalues that converge more rapidly. This is not to be expected, given the superior behavior of the Hilbert space projection method in predicting eigenvalues when δ = 0 and warrants further investigation. The Tau method clearly returns the poorest approximation to the real part of the eigenvalues and also exhibits the same stability problems, relative to convergence of the imaginary part of the eigenvalues, that were present with δ = 0. Figures (4.1)-(4.4) show convergence behavior of the imaginary and real part of the last non-zero eigenvalue for the coupled problem. Figures (4.1) and (4.2) clearly exhibit convergence and stability problems that exist when using the Tau method. Figure (4.3) shows convergence behavior with the L2 projection method with δ = 1. This projection method appeared best when δ = 1. As mentioned above, this result is not expected, and warrants further investigation. Figure (4.4) shows
N 20 40 60 80 100 120 140 160 180 200
Table 4.4: Coupled result of Tau method, δ = 0 first non-zero eigenvalue last non-zero eigenvalue 1.813141231220E-13+2.467401113167 i -2.328331971668E-12+2.986688399691E+02 i 3.381363039473E-11+2.467399221474 i -2.609234060909E-11+2.985557343955E+02 i -6.585467406376E-10+2.467488951408 i 3.122536384275E-10+2.985557371651E+02 i 5.320960028680E-10+2.466886149922 i 7.871027188168E-10+2.985557410790E+02 i 1.433304326101E-10+2.467215116948 i 1.620106678028E-09+2.985557139647E+02 i -3.238582972387E-09+2.463253463769 i 4.838114339549E-09+2.985555452318E+02 i -1.207869623498E-08+2.453631441429 i -9.425007160857E-05+2.985558906559E+02 i 1.732980771025E-08+2.483283422922 i -2.332974594969E-08+2.985559961173E+02 i 7.815160878427E-02+2.277009910907 i -1.349500070194E-03+2.985560165261E+02 i 4.517035243871E-08+2.795036974019 i -7.191520240927E-08+2.985518060386E+02 i
N 20 40 60 80 100 120 140 160 180 200
Table 4.5: Coupled result of Tau method, δ = 1 first non-zero eigenvalue last non-zero eigenvalue -3.343092260546E-02+2.467268914505 i -4.012864610018E-02+2.986687939775E+02 i -3.046995364681E-02+2.467289557063 i -3.004363359992E-02+2.985557013780E+02 i -3.046591367662E-02+2.467276642906 i -3.046422238139E-02+2.985557041466E+02 i -3.046582954389E-02+2.467768941256 i -3.046440884390E-02+2.985557039390E+02 i -3.046453325873E-02+2.468520690905 i -3.046441030609E-02+2.985556991528E+02 i -3.047390197569E-02+2.470596691565 i -3.046442146404E-02+2.985557360775E+02 i -3.890404005919E-02+2.479766768709 i -3.037677448922E-02+2.985556384643E+02 i -3.055148980721E-02+2.434383459495 i -3.046503379262E-02+2.985553356671E+02 i -2.993809503398E-02+2.270872395463 i -3.046398072003E-02+2.985569115860E+02 i 8.985638519438E-02+2.379663589928 i -3.046351060327E-02+2.985539901649E+02 i
23
Table 4.6: Coupled result of L2 inner product method based on Chebyshev basis, δ = 0 N first non-zero eigenvalue last non-zero eigenvalue 20 4.794515479078E-14+2.467401100272 i 1.264266469292E-12+2.986688399691E+02 i 40 -2.345606348042E-12+2.467401100231 i -2.268518706217E-12+2.985555697252E+02 i 60 -3.055189299250E-11+2.467401099385 i -3.028798913063E-11+2.985555702384E+02 i 80 -1.640015196130E-11+2.467401096683 i -7.921982514425E-12+2.985555702456E+02 i 100 -6.000598263179E-11+2.467401112028 i 1.471674683159E-10+2.985555702282E+02 i 120 -1.439859495773E-10+2.467401119248 i -5.120891970016E-11+2.985555702114E+02 i 140 -1.185869050967E-10+2.467401069820 i -1.209905468348E-09+2.985555702877E+02 i 160 3.639652335283E-10+2.467401247723 i -5.078418258836E-10+2.985555706296E+02 i 180 3.667867737507E-11+2.467401227815 i -9.979618652723E-10+2.985555703877E+02 i 200 3.710481723475E-10+2.467400782526 i -4.916232485453E-09+2.985555692702E+02 i
Table 4.7: Coupled result of L2 inner product method based on Chebyshev basis, δ = 1 N first non-zero eigenvalue last non-zero eigenvalue 20 -6.077228184618E-03+2.472315451479 i -7.296147769658E-03+2.986688879395E+02 i 40 -5.539132572058E-03+2.471881225063 i -5.462420024517E-03+2.985555691914E+02 i 60 -5.538374080347E-03+2.471880612819 i -5.538890964086E-03+2.985555697036E+02 i 80 -5.538374021944E-03+2.471880614097 i -5.538924742400E-03+2.985555696935E+02 i 100 -5.538373921805E-03+2.471880622592 i -5.538924666957E-03+2.985555697245E+02 i 120 -5.538374235518E-03+2.471880655294 i -5.538925170105E-03+2.985555696625E+02 i 140 -5.538373989821E-03+2.471880557778 i -5.538924641325E-03+2.985555697754E+02 i 160 -5.538373246956E-03+2.471880631905 i -5.538924178500E-03+2.985555697666E+02 i 180 -5.538374384906E-03+2.471880740717 i -5.538927030167E-03+2.985555691089E+02 i 200 -5.538375905247E-03+2.471880313327 i -5.538926914480E-03+2.985555693472E+02 i
Table 4.8: Coupled result of Hilbert Space inner product method based on Chebyshev basis, δ = 0 N first non-zero eigenvalue last non-zero eigenvalue 20 0+2.467401100186 i 0+2.991653577748E+02 i 40 0+2.467401101144 i 0+2.985555331317E+02 i 60 0+2.467401120569 i 0+2.985555331102E+02 i 80 0+2.467401172647 i 0+2.985555329554E+02 i 100 0+2.467401436480 i 0+2.985555326753E+02 i 120 0+2.467399489321 i 0+2.985555351681E+02 i 140 0+2.467403285242 i 0+2.985555305551E+02 i 160 0+2.467409597270 i 0+2.985555225387E+02 i 180 0+2.467395532838 i 0+2.985555401061E+02 i 200 0+2.467413943868 i 0+2.985555050470E+02 i
24
Tau, TwoBeams,δ = 0 −4.7211297371507754e−12+298.6691086883846i 298.68 298.66
imaginary
298.64 298.62 298.6 298.58 298.56 298.54 20
40
60
80
100
40
60
80
100
120
140
160
180
200
120
140
160
180
200
−3
1
x 10
0.5
real
0
−0.5
−1
−1.5 20
terms
Figure 4.1: Last non-zero eigenvalue, coupled problem - Tau method, δ = 0 convergence behavior with the Hilbert Space projection method with δ = 0. This projection method clearly performed best with δ = 0.
4.3
Conclusion
From the comparison of the results for the decoupled and coupled problems, we can see that the L2 inner product and Hilbert Space inner product methods are more reliable than the Tau method since they provide better accuracy and better numerical stability. This does not mean the Tau method will perform poorly for every beam equation but only for the coupled beam equations presented in this paper. The L2 inner product method gives more promising results when the damping, δ = 1 while the Hilbert Space inner product method gives better results when the damping, δ = 0. It is possible that results from the Hilbert Space inner product method would perform better when δ = 1 if we
Table 4.9: Coupled result of Hilbert Space inner product method based on Chebyshev basis, δ = 1 N first non-zero eigenvalue last non-zero eigenvalue 20 -1.398357578269E-02+2.467363269590 i -8.813979442928E-03+2.991653559500E+02 i 40 -3.113953936625E-03+2.467399282939 i -5.164412664935E-03+2.985555327638E+02 i 60 -5.650443423560E-03+2.467394987625 i -5.456026119469E-03+2.985555328427E+02 i 80 -5.539505250050E-03+2.467395494831 i -5.546996595502E-03+2.985555342605E+02 i 100 -5.538671831835E-03+2.467395605785 i -5.558298756751E-03+2.985555337636E+02 i 120 -5.538901224057E-03+2.467393663824 i -5.533184301117E-03+2.985555393982E+02 i 140 -5.537829070998E-03+2.467396294386 i -5.499261838016E-03+2.985555265588E+02 i 160 -5.536659541852E-03+2.467400659306 i -5.603798305067E-03+2.985555462832E+02 i 180 -5.530276495927E-03+2.467407454090 i -5.271875342596E-03+2.985555327163E+02 i 200 -5.529889505552E-03+2.467407865463 i -5.662180738027E-03+2.985555355425E+02 i
25
Tau, TwoBeams,δ = 1 −0.040129225959461681+298.66906269470729i 298.68 298.66
imaginary
298.64 298.62 298.6 298.58 298.56 298.54 20
40
60
80
100
40
60
80
100
120
140
160
180
200
120
140
160
180
200
−0.024 −0.026 −0.028
real
−0.03 −0.032 −0.034 −0.036 −0.038 −0.04 −0.042 20
terms
Figure 4.2: Last non-zero eigenvalue, coupled problem - Tau method, δ = 1
/UseBC, TwoBeams, δ = 1 −0.0072961477696580901+298.66888793949005i 298.68 298.66
imaginary
298.64 298.62 298.6 298.58 298.56 298.54 20
40
60
80
100
40
60
80
100
120
140
160
180
200
120
140
160
180
200
−3
−4.5
x 10
−5
real
−5.5 −6 −6.5 −7 −7.5 20
terms
Figure 4.3: Last non-zero eigenvalue, coupled problem - L2 inner product method, δ = 1
26
H6UseBC, TwoBeams, δ = 0 0+299.20938085200459i 299.3 299.2 299.1
imaginary
299 298.9 298.8 298.7 298.6 298.5 298.4 20
40
60
80
100
40
60
80
100
120
140
160
180
200
120
140
160
180
200
−6
1
x 10
0 −1
real
−2 −3 −4 −5 −6 −7 −8 20
terms
Figure 4.4: Last non-zero eigenvalue, coupled problem - Hilbert Space inner product method, δ = 0 were able to more accurately approximate the scaling parameter, ri (see equation (A.6)). In this paper, we utilized the Tau method and both L2 and Hilbert Space inner product projection methods to reduce a system of partial differential equations to a system of ordinary differential equations. The algebraic eigenvalue problem associated with the system of ordinary differential equations was then solved using standard software. The eigenvalues for the system of ordinary differential equations are approximations to eigenvalues of the original system of partial differential equations and both the L2 and Hilbert Space inner product projection methods proved to be reasonable successful in approximating these eigenvalues. It will be valuable and instructive to investigate whether these methods can be applied to higher dimensional problems.
A
Appendix
Here are some properties of Chebyshev polynomials: T0 (x) = 1 T1 (x) = x Tn (x) = 2xTn−1 (x) − Tn−2 (x), n > 2
Tn (±1) = (±1)n , n = 0, 1, . . . Qk−1 2 2 dk (n+k) m=0 (n − m ) T (±1) = (±1) , k ≥ 1, n = 0, 1, . . . n k dx (2k − 1)!!
(A.1)
Particularly, d Tn (±1) = (±1)n+1 n2 dx 2 2 d2 n+2 n (n − 1) T (±1) = (±1) n dx2 3 27
(A.2)
Define inner product Z
hu, vi = then
1 −1
√
uv dx, 1 − x2
π, hTi (x), Tj (x)i = π2 , 0,
Suppose f (x) =
N X
an Tn (x),
n=0
i=j=0 i = j 6= 0 i 6= j
(A.3)
N X dk a(k) f (x) = n Tn (x), dxk n=0
(k)
the relationship between an and an is known as recurrence rule: (k)
(k)
cn−1 an−1 = 2na(k−1) + an+1 , n ≥ 1, k ≥ 1 n
(A.4)
where cj =
(
2, 1
j=0 j = 1, 2, . . .
PN (1) If we suppose aN = i=0 sn,i ai , n = 0, . . . , N then by substituting this sum into A.4 and equating coefficients of ai , it is fairly straightforward to show that the sn,i have the following properties: sn,i sn−1,n sn−1,i
= 0 , 2n = , cn−1 sn+1,i , = cn−1
n = 0, . . . , N, i = 0, . . . , n n = N, . . . , 1 n = N − 1, . . . , 1, i = n + 2, . . . , N .
Thus if we let M = [sn,i ] then from A.4 and mathematical induction, it can be shown that (k) a0 a0 (k) a1 a1 . = Mk .. . . . . (k) aN a
(A.5)
N
In order to normalize the maximum coefficient of u1 ,v1 ,u2 ,v2 in 3.19 to be approximately 1, we introduce the scale parameter γn through the following normalization, Z 1 2 (G′′n (x)) 1 √ dx = 1 γn2 −1 1 − x2 where Gn (x) is defined in equation (3.7). It is not possible to get a closed-form expression for γn . Therefore, we obtained approximate values to γn by performing a discrete least squares fit to the data generated from ri = a(1 + i)b = γi , i = 0, . . . , 25 . In this way we find that a ≈ 4.98479 and b ≈ 2.00304 and thus that ri = 4.98479(1 + i)2.00304 . which is still a good approximation to γi for i ≈ 200. 28
(A.6)
B B.1
Matlab Code main Listing B.1: main.m
1
% main f i l e
2 3 4 5
clear ; clear a l l ; close a l l ;
6 7 8 9
% N i s b e t w e e n n1 and n2 n1 = 2 0 ; n2 = 2 0 0 ;
10 11 12 13
% keep only e i g e n v a l u e s of l e n g t h < radius with % non−n e g a t i v e i m a g i n a r y p a r t radius = 300;
14 15 16 17 18
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
% %% %% Tau method : Decoupled problem k ( x ) =0, 1 beam f o r p = 1 %boundary c o n d i t i o n : 1 or 2 b u t i n t h i s problem 1 & 2 a r e same conditions f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = Tau1Beam ( n1 , r , p ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = Tau1Beam (m, r , p ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave1 ( 0 , r e s ( : ,m) , r , p , 0 , n1 , n2 ) ; end end end
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
%% Tau method : Coupled problem k ( x ) n o t e q u a l s t o 0 so t h a t 2 beams f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = Tau2Beam ( n1 , r ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = Tau2Beam (m, r ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave2 ( 0 , r e s ( : ,m) , r , 1 , n1 , n2 ) ; end end
58 59 60
%% L2 i n n e r p r o d u c t b a s e d on Chebshev b a s i s : %% v1 & v2 u s e boundary c o n d i t i o n : Decoupled problem (1 beam )
29
61
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
f o r p = 1 %boundary c o n d i t i o n : 1 or 2 ; b u t i n t h i s problem 1 & 2 a r e same conditions f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = L2UseBC1Beam ( n1 , r , p ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = L2UseBC1Beam (m, r , p ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave1 ( 1 , r e s ( : ,m) , r , p , 0 , n1 , n2 ) ; end end end
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
%% L2 i n n e r p r o d u c t b a s e d on Chebshev b a s i s : %% v1 & v2 u s e boundary c o n d i t i o n : c o u p l e d problem (2 beams ) f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = L2UseBC2Beam ( n1 , r ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = L2UseBC2Beam (m, r ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave2 ( 1 , r e s ( : ,m) , r , 1 , n1 , n2 ) ; end end
102 103 104 105 106
107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
%% H e l b e r t Space I n n e r Product b a s e d on Chebshev b a s i s : Decoupled problem f o r p = 1 %boundary c o n d i t i o n : 1 or 2 ; b u t i n t h i s problem 1 & 2 a r e same conditions f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = HSUseBC1Beam ( n1 , r , p ) ;%% H Space tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = HSUseBC1Beam (m, r , p ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave1 ( 2 , r e s ( : ,m) , r , p , 0 , n1 , n2 ) ; end end end
126
30
127
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
B.2
%% H e l b e r t Space I n n e r Product b a s e d on Chebshev b a s i s : Decoupled problem = 2 beams f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = HSUseBC2Beam ( n1 , r ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = HSUseBC2Beam (m, r ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave2 ( 2 , r e s ( : ,m) , r , 1 , n1 , n2 ) ; end end
c Listing B.2: c.m
1 2 3 4 5 6 7 8
B.3
% c(k) function r = c ( k ) i f ( k==0) r = 2.0; else r = 1.0; end end
ChebyshevExpansion Listing B.3: ChebyshevExpansion.m
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
B.4
% compute t h e c o e f f i c i e n t s o f t h e Chebyshev e x p a n s i o n s % f o r the given func function r = ChebyshevExpansion ( func , n ) r=zeros ( 1 , n+1) ; f1 = func (1) ; f 2 = f u n c ( −1) ; t =1:n−1; f o r m=0:n tmp = ( f 1 + f 2 ) / 2 ; tmp = tmp + f u n c ( cos ( pi ∗ t /n ) ) ∗ cos ( pi ∗m∗ t /n ) ’ ; i f (m==0 | | m==n ) r (m+1) = tmp/n ; else r (m+1) = 2∗tmp/n ; end f 2 = −f 2 ; end end
ChebyshevMultiplication Listing B.4: ChebyshevMultiplication.m
1 2 3 4 5 6
% f o r m u l a f o r f ( x ) g ( x ) where f ( x ) and g ( x ) a r e e x p a n s i o n s % i n term o f Chebyshev b a s i s function r = C h e b y s h e v M u l t i p l i c a t i o n ( c1 ) L = length ( c1 ) −1; r = zeros ( 2 ∗L+1 ,L+1) ; for n = 0 : L
31
for m = 0 : L r (m+n+1 , m+1) = r (m+n+1 , m+1) + c1 ( n+1) / 2 ; i f m>=n k = m − n; else k = n − m; end r ( k+1 , m+1) = r ( k+1 , m+1) + c1 ( n+1) / 2 ; end
7 8 9 10 11 12 13 14 15
end
16 17
end
B.5
GetClosestPoint Listing B.5: GetClosestPoint.m
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
% c h o o s e a complex number from ptmp which i s c l o s e s t t o a p o i n t function r = G e t C l o s e s t P o i n t ( ptmp , a p o i n t ) t = 0.01; l = 0; while l
A Project Submitted to the Faculty of the Department of Mathematics and Statistics University of Minnesota, Duluth by
Zhengfei Rui In partial fulfillment of the requirements for the degree of Master of Science in Applied and Computational Mathematics
Advisor: Steven A. Trogdon
June 24, 2014
Acknowledgments First and foremost, I would like to show my deepest gratitude to my supervisor, Dr. Steven A. Trogdon, a respectable, responsible and resourceful professor, who has provided me with valuable guidance in every stage of the writing of this project. Without his enlightening instruction, impressive kindness and patience, I could not have completed it. His keen and vigorous academic observation enlightens me not only in this paper work but also in my future study. Secondly, I shall extend my thanks to Dr. Zhuangyi Liu, Dr. Jiann-Shiou Yang and Dr. Debao Zhou for all their kindness and help. I would also like to thank all my teachers who have helped me to develop the fundamental and essential academic competence. Last but not least, I’ d like to thank my elder brother Zhenhua Rui for his encouragement and support. Zhengfei Rui
1
Abstract Chebyshev polynomials may be used to construct basis functions which satisfy fairly general boundary conditions. For the purpose of this report, these boundary conditions are associated with finding solutions to a system of coupled partial differential equations in space and time. These basis functions along with L2 and Hilbert space inner products may be used to spatially discretize the system of partial differential equations yielding a system of ordinary differential equations in time. By exploiting properties of Chebyshev polynomials the system of ordinary differential equations may be converted to an algebraic eigenvalue problem which may be solved using standard software. Using basis functions that satisfy boundary conditions is shown to give more accurate results for eigenvalues than when the Tau method, where basis function only approximately satisfy boundary conditions, is used.
2
Contents Acknowledgments
1
Abstract
2
List of Tables
4
List of Figures
5
List of Matlab Code
6
1 Introduction
7
2 Notation and Formulae
8
3 Beam Equations 3.1 Construction of a Chebyshev Basis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 Spatially Discretized Equations Using a L2 Inner Product . . . . . . . . . . . . . . . 3.3 Spatially Discretized Equations Using a Hilbert Space Inner Product . . . . . . . . .
10 10 11 14
4 Numerical results 4.1 Decoupled problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Coupled problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
21 21 21 25
A Appendix
27
B Matlab Code B.1 main . . . . . . . . . . . B.2 c . . . . . . . . . . . . . B.3 ChebyshevExpansion . . B.4 ChebyshevMultiplication B.5 GetClosestPoint . . . . . B.6 DrawSave1 . . . . . . . B.7 DrawSave2 . . . . . . . B.8 Tau1Beam . . . . . . . . B.9 Tau2Beam . . . . . . . . B.10 L2NB1Beam . . . . . . B.11 L2NB2Beam . . . . . . B.12 HSNB1Beam . . . . . . B.13 HSNB2Beam . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
3
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
29 29 31 31 31 32 32 33 34 35 36 37 38 40
List of Tables 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9
Decoupled result of Tau method . . . . . . . . . . . . . . . . . . . . . . . . . . . . Decoupled result of L2 inner product method based on Chebyshev basis . . . . . . Decoupled result of Hilbert Space inner product method based on Chebyshev basis Coupled result of Tau method, δ = 0 . . . . . . . . . . . . . . . . . . . . . . . . . Coupled result of Tau method, δ = 1 . . . . . . . . . . . . . . . . . . . . . . . . . Coupled result of L2 inner product method based on Chebyshev basis, δ = 0 . . . Coupled result of L2 inner product method based on Chebyshev basis, δ = 1 . . . Coupled result of Hilbert Space inner product method based on Chebyshev basis, δ Coupled result of Hilbert Space inner product method based on Chebyshev basis, δ
4
. .
22 22 22 . 23 . 23 . 24 . 24 = 0 24 = 1 25
List of Figures 4.1 4.2 4.3 4.4
Last Last Last Last
non-zero non-zero non-zero non-zero
eigenvalue, eigenvalue, eigenvalue, eigenvalue,
coupled coupled coupled coupled
problem problem problem problem
5
-
Tau method, δ = 0 . . . . . . . . . . Tau method, δ = 1 . . . . . . . . . L2 inner product method, δ = 1 . . Hilbert Space inner product method,
. . . δ
. 25 . 26 . 26 = 0 27
List of Matlab Code B.1 B.2 B.3 B.4 B.5 B.6 B.7 B.8 B.9 B.10 B.11 B.12 B.13
main.m . . . . . . . . . . . c.m . . . . . . . . . . . . . . ChebyshevExpansion.m . . ChebyshevMultiplication.m GetClosestPoint.m . . . . . DrawSave1.m . . . . . . . . DrawSave2.m . . . . . . . . Tau1Beam.m . . . . . . . . Tau2Beam.m . . . . . . . . L2NB1Beam.m . . . . . . . L2NB2Beam.m . . . . . . . HSNB1Beam.m . . . . . . . HSNB2Beam.m . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
6
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
. . . . . . . . . . . . .
29 31 31 31 32 32 33 34 35 36 37 38 40
1
Introduction
The spectral method is widely used to solve partial differential equations (PDEs). The method consists in approximating the solution to a PDE as a linear combination of a set of orthogonal basis functions. This approximation is called a spectral approximation. Once the approximation is substituted into the PDE, a projection procedure is employed to convert the original PDE into a system of equations for coefficients that multiply the basis functions in the spectral approximation. Chebyshev polynomials are often chosen as the basis functions used in implementing the spectral method. These polynomials are defined as Tn (x) = cos(n arccos x), −1 ≤ x ≤ 1, n = 0, 1, . . . and are orthogonal with respect to the inner product Z 1 uv √ dx hu, vi = 1 − x2 −1
(1.1)
(1.2)
In general, for a PDE of order 2k Lu = g, −1 ≤ x ≤ 1
(1.3)
with boundary conditions Bl,i u(−1) = 0, Br,i u(1) = 0, i = 1, . . . k
(1.4)
where L, Bl,i , Br,i are linear differential operators and g is a known function of x, we look for an approximate solution in the N + 1 dimensional Chebyshev function space VN = span{Tn (x), n = 0, . . . N } as u ∼ uN =
N X
an Tn (x)
(1.5)
(1.6)
n=0
and require that LuN − g be orthogonal to VN with respect to our inner product. This orthogonality or projection requirement means that hLuN − g, Ti (x)i = 0, i = 0, . . . N
(1.7)
which is equivalent to N X
n=0
an hLTn (x), Ti (x)i = hg, Ti (x)i, i = 0, . . . N
(1.8)
To satisfy the boundary conditions (1.4), we require that Bl,i uN (−1) = 0, Br,i uN (1) = 0, i = 1, . . . k
(1.9)
which is equivalent to N X
n=0
an Bl,i Tn (−1) = 0,
N X
an Br,i Tn (1) = 0, i = 1, . . . k
(1.10)
n=0
In equations (1.8) and (1.10), there are N + 1 + 2k equations with N + 1 unknowns. To make the number of equations and unknowns consistent, we apply the inner product in (1.8) for i = 0, . . . , N − 2k; the remaining 2k equations in (1.8) are dropped from the analysis.
7
2
Notation and Formulae For convenience, we introduce some notation. Suppose
f (x) =
N X
fi Ti (x) = f T TN
(2.1)
i=0
where
T0 (x) f0 f = ... , TN = ...
(2.2)
TN (x)
fN
and f T denotes the transpose of f . From now on, f (x) = f T TN is denoted as
V
N f f (x) −→
(2.3)
where VN is the Chebyshev function space defined in the Introduction (1.5). This notion means that f (x) may be expanded in the function space VN as the vector f . Furthermore, using properties of Chebyshev polynomials, it can be shown that V
N f⇒ f (x) −→
dk VN M kf , k ≥ 1 f (x) −→ k dx
(2.4)
where the matrix M is of order (N + 1) × (N + 1) and is defined in the Appendix (A.5). Now let g0 VN .. g(x) −→ g = . .
(2.5)
gN
By using the multiplication law,
Ti (x)Tj (x) =
1 1 Ti+j (x) + T|i−j| (x) 2 2
(2.6)
for the Chebyshev polynomials, the product f (x)g(x) may be approximated as f (x)g(x)
= = = =
N X
fi Ti (x)
i=0 N N X X
i=0 j=0 N N X X
N X
fi gj Ti (x)Tj (x) fi g j
i=0 j=0 2N X
1 2
gj Tj (x)
j=0
1 1 Ti+j (x) + T|i−j| (x) 2 2
S1,i Ti (x) +
i=0
N
1X S2,i Ti (x) 2 i=0
where S1,i =
X
fp g p ,
p+q =i 0 ≤ p, q ≤ N
8
0 ≤ i ≤ 2N
(2.7)
and X
S2,i =
0≤i≤N
fp g p ,
|p − q| = i 0 6 p, q 6 N
(2.8)
From (2.7) we have that S1,0 S1,1 S1,2 .. .
= = =
f0 g 0 f0 g 1 + f1 g 0 f0 g 2 + f1 g 1 + f2 g 0
S1,N
=
f0 gN + f1 gN −1 + · · · + fN g0
and thus
S1,0 f0 S1,1 f1 .. = .. . . S1,N
f0
fN
fN −1
···
Likewise from (2.8) we have that S2,0 S2,1 .. . S2,N
= f0 g 0 + f1 g 1 + · · · + fN g N = f0 g1 + f1 g2 + · · · + fN −1 gN + f1 g0 + f2 g1 + · · · + fN gN −1
(2.9)
(2.10) (2.11) (2.12) (2.13)
= f0 g N + fN g 0
and thus that S2,0 f0 S2,1 f1 .. = .. . . S2,N
g 0 .. . ≡ W1 g. .. . f0 gN
fN
f1 f0 + f2
f2 f1 + f3
··· ···
fN −1 fN −2 + fN
0
···
···
0
Hence we get that
f (x)g(x)
= = = =
g 0 .. fN −1 . ≡ W2 g. .. . f0 gN fN
(2.14)
2N N 1X 1X S1,i Ti (x) + S2,i Ti (x) 2 i=0 2 i=0 N N 2N 1X 1 X 1X S1,i Ti (x) + S2,i Ti (x) + S1,i Ti (x) 2 i=0 2 i=0 2 i=N +1 1 1 1 T T T [W1 g] TN + [W2 g] TN + [W3 g] TN∗ 2 2 2 1 1 T T [(W1 + W2 )g] TN + [W3 g] TN∗ 2 2
where TN∗
TN +1 (x) .. = . .
(2.15)
T2N (x)
If we let Wf = 21 (W1 + W2 ), and Wf∗ = 12 W3 , then we get that T T f (x)g(x) = [Wf g] TN + Wf∗ g TN∗
(2.16)
where Wf∗ is not given explicitly because it consists of only terms of order higher than N , which become 0 after taking the inner product with Tk (x), k = 0, 1, . . . N. 9
3
Beam Equations
We are given the following system of partial differential equations describing the motion of two coupled beams ∂ 4 u1 ∂u1 ∂ 2 u1 = −p1 − k(x) (u1 − u2 ) − δk(x) , 2 ∂t ∂x4 ∂t ∂ 2 u2 ∂ 4 u2 ∂u2 = −p + k(x) (u1 − u2 ) − δk(x) . 2 ∂t2 ∂x4 ∂t
(3.1a) (3.1b)
with boundary conditions
2
u1 (t, −1) =u1 (t, 1) = 0,
(3.2a)
2
∂ u1 ∂ u1 (t, −1) = (t, 1) = 0, 2 ∂x ∂x2 u2 (t, −1) =u2 (t, 1) = 0,
(3.2b) (3.2c)
2
2
∂ u2 ∂ u2 (t, −1) = (t, 1) = 0. ∂x2 ∂x2
(3.2d)
where p1 > 0, p2 > 0, and δ ≥ 0 are given parameters and k(x) is a known function. ∂u2 1 If we introduce the substitutions, v1 = ∂u ∂t and v2 = ∂t , then the original equations (3.1a) and (3.1b) become ∂u1 ∂t ∂v1 ∂t ∂u2 ∂t ∂v2 ∂t
(3.3a)
= v1 , 4
= −p1
∂ u1 − k(x) (u1 − u2 ) − δk(x)v1 , ∂x4
(3.3c)
= v2 , = −p2
(3.3b)
∂ 4 u2 + k(x) (u1 − u2 ) − δk(x)v2 . ∂x4
(3.3d)
Next, in section 3.1 we will construct basis functions that satisfy the boundary conditions using Chebyshev polynomials. Then, in sections 3.2 and 3.3 we use two different projection methods to set up eigenvalue (linear algebra) problems which can be solved, and in section 4 we compare the numerical results these two different projection methods give us.
3.1
Construction of a Chebyshev Basis
In this section we construct a set of functions that span a vector space GN , where each basis element of GN satisfies boundary conditions (3.2a)-(3.2d). We then expand u1 ,v1 ,u2 ,v2 in GN as they should satisfy the associated boundary conditions, and solve (3.3a)-(3.3d) with boundary conditions (3.2a)-(3.2d). Let GN = span{Gn (x), n = 0, ...N } where Gn = Tn (x) + a1 Tn+1 (x) + a2 Tn+2 (x) + a3 Tn+3 (x) + a4 Tn+4 (x) = Tn (x) +
4 X
ai Tn+i (x)
i=1
10
Plugging the above Gn (x) into (3.2a) and (3.2b), and using formula (A.2) from Appendix A, we get
Gn (−1) = 0 ⇒ 1 + Gn (1) = 0 ⇒ 1 +
4 X
(−1)i ai = 0,
(3.5a)
ai = 0,
(3.5b)
i=1
4 X i=1
4
2 2 d2 n2 (n2 − 1) X i (n + i) [(n + i) − 1]ai G (−1) = 0 ⇒ (−1) + = 0, n dx2 3 3 i=1
(3.5c)
4
n2 (n2 − 1) X (n + i)2 [(n + i)2 − 1]ai d2 (−1)i Gn (−1) = 0 ⇒ + = 0. 2 dx 3 3 i=1
(3.5d)
These are a set of algebraic equations with four unknowns, ai , i = 1, ...4. We use Mathematica to solve for these unknowns and get that a1 (n) = a3 (n) = 0, a2 (n) = −
2(n + 2)(2n2 + 8n + 15) (n + 1)(2n2 + 4n + 3) , a (n) = . 4 (n + 3)(2n2 + 12n + 19) (n + 3)(2n2 + 12n + 19)
(3.6)
Therefore, Gn (x) = Tn (x) + a2 TN +2 (x) + a4 Tn+4 (x).
(3.7)
It can be shown that G0 (x), G1 (x), ...GN (x) are linear independent, and thus span the subspace GN of VN +4 . It is easy to see that, every g(x) ∈ GN is representable in VN +4 . In fact, we have that VN +4
G
N g ⇒ g(x) −→ Gg g(x) −→
(3.8)
where
1
0 a2 (0) 0 a4 (0) G= 0 . .. . .. . .. 0
0
0
1
0
0
1
a2 (1) 0 a4 (1)
0 a2 (2) 0
.. . .. . 0
a4 (2) 0 .. . 0
···
0 .. . .. . .. .
0 . 1 0 a2 (N ) 0
(3.9)
a4 (N )
Thus G has order (N + 5) × (N + 1) and represents the linear transformation from GN to VN +4 .
3.2
Spatially Discretized Equations Using a L2 Inner Product
Let’s map the u1 , u2 , v1 , v2 to the VN +4 , so that we can apply the product formula (2.16) and inner product relation (A.3) to them with Tn (x), n = 0, . . . , N . We also expand k(x) in VN +4 since it will
11
multiply functions from VN +4 instead of GN . The expansions are as follows: V N +4 a0 (t) u1 −→ Ga(t) GN VN +4 1 a(t) = ... ⇒ ∂u u1 −→ ˙ ∂t −→ Ga(t) ∂ 4 u1 VN +4 4 aN (t) ∂x4 −→ M Ga(t), VN +4 b0 (t), u2 −→ Gb(t) GN VN +4 2 ˙ b(t) = ... ⇒ ∂u u2 −→ ∂t −→ Gb(t) 4 ∂ u2 VN +4 4 bN (t) ∂x4 −→ M Gb(t) c0 (t), N +4 v V−→ Gc(t) GN 1 c(t) = ... ⇒ v1 −→ N +4 ∂v1 V−→ Gc(t) ˙ ∂t cN (t) d0 (t) N +4 v V−→ Gd(t) GN .. 2 v2 −→ d(t) = . ⇒ N +4 ∂v2 V−→ ˙ Gd(t), ∂t dN (t) k0 . VN +4 k(x) −→ k = ... , kN +4
where the superposed dot denotes time differentiation. Upon substituting from the above expansions into (3.3a), we get that N X
a˙i (t)Gi (x) =
N X
ci (t)Gi (x)
i=0
i=0
and since the Gi (x), i = 0, . . . , N are linearly independent, we deduce that a˙ i (t) = ci (t), i = 0, . . . , N. This can be written in matrix-vector form as IN +1 a(t) ˙ = IN +1 c(t)
(3.11)
where IN +1 is the identity matrix of order N +1. This is the spatially discretized equation associated with (3.3a). Likewise, we can also get the spatially discretized equation associated with (3.3c) as ˙ = IN +1 d(t). IN +1 b(t)
(3.12)
The remaining equations to be discretized are equations (3.3b) and (3.3d). Using (2.3) for functions expanded in the function space, VN +4 we get that k(x)u1 (x) = [Wk Ga(t)]T TN +4 + [Wk∗ Ga(t)]T TN∗ +4 , k(x)u2 (x) = [Wk Gb(t)]T TN +4 + [Wk∗ Gb(t)]T TN∗ +4 , k(x)v1 (x) = [Wk Gc(t)]T TN +4 + [Wk∗ Gc(t)]T TN∗ +4 , k(x)v2 (x) = [Wk Gd(t)]T TN +4 + [Wk∗ Gd(t)]T TN∗ +4 . With the function space expansion of u1 , u2 and v1 in VN +4 , we then obtain from (3.3b) that [Gc(t)] ˙ T TN +4 = [−p1 M 4 Ga(t)]T TN +4 − [Wk Ga(t)]T TN +4 − [Wk∗ Ga(t)]T TN∗ +4 ,
+[Wk Gb(t)]T TN +4 + [Wk∗ Gb(t)]T TN∗ +4 − δ[Wk Gc(t)]T TN +4 − δ[Wk∗ Gc(t)]T TN∗ +4 . 12
After taking the inner product of these equations with Tj (x), j = 0, . . . N + 4 we are left with N + 5 equations in N + 1 unknowns. To make the system of equations consistent, we keep the first N + 1 equations and drop the last 4 equations. In fact, if k(x) = 0, i.e. the two beams are independent of each other, the right hand side of the last 4 equations is identically zero which makes the system of equation inconsistent unless those equations are dropped. Upon taking the inner product with Tj (x), j = 0, . . . N , we get ˙ = [−p1 M 4 G]j a(t) − [Wk G]j a(t) + [Wk G]j b(t) − δ[Wk G]j c(t), j = 0, . . . N. [G]j c(t) This can be written in matrix-vector form as ˙ = A1 a(t) − Ba(t) + Bb(t) − δBc(t). G1 c(t)
(3.13)
Equations (3.13) are the spatially discretized equations corresponding to (3.3b), where G1 , of order(N + 1) × (N + 1) is the submatrix obtained from the first N + 1 rows of G. And A1 , of order(N + 1) × (N + 1) is the submatrix obtained from the first N + 1 rows of −p1 M 4 G. Also B, of order(N + 1) × (N + 1) is the submatrix obtained from the first N + 1 rows of Wk G. In the same way, we get ˙ = [−p2 M 4 G]j b(t) + [Wk G]j a(t) − [Wk G]j b(t) − δ[Wk G]j d(t), j = 0, . . . N [G]j d(t) This can also be written in matrix-vector form as, ˙ = A2 b(t) + Ba(t) − Bb(t) − δBd(t). G1 d(t)
(3.14)
where A2 , of order(N +1)×(N +1) is the submatrix obtained from the first N +1 rows of −p2 M 4 G. If we let a(t) c(t) (3.15) z(t) = b(t) d(t)
then (3.11)-(3.14) may be combined as the following
IN +1 0 0 0
0 0 G1 0
0 IN +1 0 0
0 0 0 0 z(t) ˙ = A1 − B 0 B G1
IN +1 0 −δB 0
0 0 B A2 − B
which may more compactly be written as
˙ Qz(t) = Dz(t).
0
IN +1 z(t) 0 −δB (3.16)
The above equations are the spatially discretized system of equations associated with using basis functions which satisfying the boundary conditions. The eigenvalue problem for these equations is ˆ where obtained by setting z(t) = eλt z, ˆ a ˆb0 a ˆ0 cˆ0 dˆ0 cˆ . ˆ = .. , cˆ = .. , dˆ = .. . ˆ= zˆ = (3.17) .. , b . . . ˆ , a b ˆ ˆ a ˆN cˆN bN dN dˆ Therefore we obtain
λQeλt zˆ = Deλt zˆ ⇒ λzˆ = Q−1 D zˆ which is the eigenvalue problem corresponding to the equations (3.16). 13
(3.18)
3.3
Spatially Discretized Equations Using a Hilbert Space Inner Product
The Hilbert space inner product that we use in this section requires that we calculate the L2 inner product, hD 2 Gi , D 2 Gj iL2 . From equations (2.4) and (A.5) it is fairly straightforward to show that D 2 Gi becomes quite large as the index, i becomes large and thus elements in the matrix representing the above inner product can become quite large. It is the u1 and u2 in equations (3.3a)(3.3d) that are needed in calculating the elements of hD 2 Gi , D 2 Gj iL2 . Thus, in order to insure that all variables have the same relative order of magnitude, we introduce a scaling parameter, ri to normalize the coefficients of the vectors and let ai (t) ri u1 v 1 GN ci (t) (3.19) z= , u2 −→ bi (t) v2 ri di (t) where the following series expansions
u1 =
N X
ai (t)
i=0
v1 =
N X
Gi (x) , ri
ci (t)Gi (x),
i=0
u2 =
N X
bi (t)
i=0
v2 =
N X
Gi (x) , ri
di (t)Gi (x)
i=0
are used and where ri , the scaling parameter is defined in equation (A.6) in Appendix A. Upon substituting the above expressions for u1 , v1 , u2 , v2 into (3.3a)-(3.3d), we get that N X
N
Gi (x) X a˙i (t) ci (t)Gi (x), = ri i=0 i=0 N X i=0
N X
c˙i (t)Gi (x) = −p1
Gi (x) = b˙i (t) ri i=0 N X i=0
N X
N X
ai (t)
i=0
(3.20)
N N X X Gi (x) D 4 Gi (x) (ai (t) − bi (t)) ci (t)Gi (x), − k(x) − δk(x) ri ri i=0 i=0
(3.21)
(3.22)
di (t)Gi (x),
i=0
d˙i (t)Gi (x) = −p2
N X i=0
bi (t)
N N X X D 4 Gi (x) Gi (x) (ai (t) − bi (t)) di (t)Gi (x). + k(x) − δk(x) ri ri i=0 i=0
(3.23)
14
Now let the left hand side (LHS) of the above equations be N X
Gi (x) a˙ i (t) ri i=0 N X c˙i (t)Gi (x) i=0 . zl = N X G (x) i b˙ i (t) r i i=0 X N d˙i (t)Gi (x)
(3.24)
i=0
from which we may define the following four vectors as a basis for representing the elements of zl
zˆ(1)
Gj (x) 0 0 0 rj 0 (2) Gj (x) (3) (4) 0 = 0 , zˆ = 0 , zˆ = Gj (x) , zˆ = 0 . 0 rj 0 Gj (x) 0 0
(3.25)
In order to obtain the best approximation, we will utilize a suitably defined Hilbert space inner product [1, 4] to spatially discretize the system of equations given above. We provide the following theorems that relate to the system of equations in section 3 as justification for the use of a Hilbert space inner product: Theorem 3.1 The operator A generates a semigroup, eAt , on the Hilbert space H. Therefore, the system is well-posed, i.e, for any Z0 ∈ H, the unique solution to the system is Z(t) = eAt Z0 , t > 0. [4] Theorem 3.2 The semigroup eAt is exponentially stable. [4] where H = U1 × V 1 × U2 × V 2 , U1 = H01 (−1, 1) ∩ H 2 (−1, 1) , V1 = L2 (−1, 1),
and
U2 = H01 (−1, 1) ∩ H 2 (−1, 1) , V2 = L2 (−1, 1), 1
f 2 (x) √ L (−1, 1) = f (x) | hf, f iL2 ≡ dx < ∞ , 1 − x2 −1 n o H k (−1, 1) = f (x) | f (x), . . . , f (k−1) (x) ∈ L2 (−1, 1) , H01 (−1, 1) = f (x) ∈ H 1 (−1, 1) | f (−1) = f (1) = 0 . 2
Z
Thus, relative the beam equations in section 3, for any z = (u1 , v1 , u2 , v2 )T ∈ H and any ˆ z = (ˆ u1 , vˆ1 , u ˆ2 , vˆ2 )T ∈ H, we may define the following inner product:
ˆ H = (u1 , v1 , u2 , v2 )T , (ˆ hz, zi u1 , vˆ1 , u ˆ2 , vˆ2 )T H ≡
p1 hD 2 u1 , D 2 u ˆ1 iL2 + p2 hD 2 u2 , D 2 u ˆ2 iL2 + hv1 , vˆ1 iL2 + hv2 , vˆ2 iL2 +hk(x)(u1 − u2 ), (uˆ1 − u ˆ2 )iL2 .
15
We now compute hzl , zˆ(1) iH , hzl , zˆ(2) iH , hzl , zˆ(3) iH and hzl , zˆ(4) iH , where the simplified notation, h·, ·i ≡ h·, ·iL2 is employed: N N a˙ i − b˙i X X a ˙ i hD 2 Gi , D 2 Gj i + hk(x)Gi , Gj i, (3.26) hzl , zˆ(1) iH = p1 r r ri rj i j i=0 i=0 hzl , zˆ(2) iH = hzl , zˆ(3) iH
N X i=0
(3.27)
c˙i hGi , Gj i,
N N a˙ i − b˙i X X ˙bi = p2 hD 2 Gi , D 2 Gj i − hk(x)Gi , Gj i, r r ri rj i j i=0 i=0
hzl , zˆ(4) iH =
N X i=0
d˙i hGi , Gj i.
(3.28)
(3.29)
Now let
N X
ci Gi (x) i=0 N N N 4 X X X D Gi (x) Gi (x) ai (ai − bi ) ci Gi (x) − k(x) − δk(x) −p1 r r i i i=0 i=0 i=0 zr = N X di Gi (x) i=0 N N N X X X D 4 Gi (x) Gi (x) −p2 (ai − bi ) di Gi (x) bi + k(x) − δk(x) ri ri i=0 i=0 i=0
(3.30)
be the right hand side (RHS) of our system of equation and compute hzr , zˆ(1) iH , hzr , zˆ(2) iH , hzr , zˆ(3) iH and hzr , zˆ(4) iH , N N c˙i − d˙i X X c ˙ i hzr , zˆ(1) iH = p1 hD 2 Gi , D 2 Gj i + hGi , Gj i, (3.31) r rj j i=0 i=0 ˙ N N N a ˙ − b X X X i i a˙ i hzr , zˆ(2) iH = −p1 ci hk(x)Gi , Gj i, (3.32) hD 4 Gi , Gj i − hk(x)Gi , Gj i − δ r ri i=0 i i=0 i=0 ˙i N ˙ N c ˙ − d X X i di hzr , zˆ(3) iH = p2 hD 2 Gi , D 2 Gj i − hG, Gj i, (3.33) r rj j i=0 i=0 N N N ˙ a˙ i − b˙i X X X b i di hk(x)Gi , Gj i. (3.34) hD 4 Gi , Gj i + hk(x)Gi , Gj i − δ hzr , zˆ(4) iH = −p2 r ri i=0 i=0 i=0 i
16
The inner product operation applied to the LHS should be equal to that applied to the RHS, thus ˙ N N a˙ i (t) − bi (t) X X a˙ i (t) 2 p1 hD Gi , D 2 Gj i + hk(x)Gi , Gj i r r ri rj i j i=0 i=0 ˙i (t) N N c ˙ (t) − d X X i c˙i (t) 2 hD Gi , D 2 Gj i + hGi , Gj i, (3.35a) = p1 r rj j i=0 i=0 ˙i (t) N N N a ˙ (t) − b X X X i a˙ i (t) 4 c˙i (t)hGi , Gj i = −p1 hD Gi , Gj i − hk(x)Gi , Gj i ri ri i=0 i=0 i=0 N X
−δ p2
N ˙ X bi (t) i=0
ri rj
hD 2 Gi , D 2 Gj i − = p2
i=0
d˙i (t)hGi , Gj i = −p2
ri rj
i=0
N ˙ X di (t)
rj
i=0
N X
N a˙ i (t) − b˙i (t) X
ri
hD 4 Gi , Gj i +
(3.35b)
ci (t)hk(x)Gi , Gj i,
hk(x)Gi , Gj i
hD 2 Gi , D 2 Gj i −
N ˙ X bi (t) i=0
i=0
N c˙i (t) − d˙i (t) X rj
i=0
N a˙ i (t) − b˙i (t) X ri
i=0
−δ
N X i=0
(3.35c)
hGi , Gj i,
hk(x)Gi , Gj i (3.35d)
di (t)hk(x)Gi , Gj i.
Here we will perform some calculations to simplify items in the above equations. For every g(x) ∈ GN , g(x) can expressed in terms of either Gi or Ti as follows: g(x) =
N X
gi Gi (x) =
N +4 X
(3.36)
gˆi Ti (x), i = 0, . . . , N + 4.
i=0
i=0
From (3.8), we know that gˆ = Gg so that gˆi =
N X
(3.37)
Gi,k gk , i = 0. . . . , N + 4,
k=0
and g(x) =
N N +4 X X
Gi,k gk Ti (x) =
i=0 k=0
N X
k=0
N +4 X
!
(Gi,k Ti (x)) gk =
i=0
N X i=0
N +4 X
!
(Gk,i Tk (x)) gi .
k=0
(3.38)
Comparing the above to (3.36), it is easy to see that Gi (x) =
N +4 X
(3.39)
Gk,i Tk (x), i = 0, . . . , N.
k=0
So if we define GN = GT TN +4 , where G is (N + 5) × (N + 1) and TN +4 is (N + 5) × 1, then from (3.36) g(x) = g T GN . Likewise, for the pth derivative g
(p)
(x) =
N X i=0
(p) gi Gi (x)
=
N +4 X i=0
17
(p) gˆi Ti (x)
=
N +4 X i=0
(p)
gˆi Ti (x)
(3.40)
where from equation (2.4) (p)
gˆi
=
N +4 X
(M p )i,k gˆk .
(3.41)
k=0
So we can express g (p) (x) as follows: +4 N +4 N X X
g (p) (x) =
(M p )i,k gˆk Ti (x)
i=0 k=0
+4 N +4 N X X
=
(M p )i,k
=
+4 N +4 N X X
m=0 N X
=
Gk,m gm Ti (x)
m=0
i=0 k=0
N X
N X
!
(M p )i,k Gk,m Ti (x) gm
i=0 k=0
G(p) m (x)gm .
m=0
Now if we switch the indices m and i, we get that g (p) (x)
=
N P
i=0 (p)
Gi (x) =
+4 NP +4 NP m=0 k=0
(p)
Gi (x)gi , where
(p)
(M p )m,k Gk,m Tm (x), i = 0, ..., N . Therefore if we define GN = (M p G)T TN +4
we may write (3.40) as (p)
g (p) (x) = g T GN .
(3.42)
Now we simplify the other inner products that appear in equations (3.35a) - (3.35d). *N +4 + N +4 N +4 N +4 X X X X Gk,i Gm,j hTk (x), Tm (x)i hGi , Gj i = Gm,j Tm (x) = Gk,i Tk (x), k=0 m=0
m=0
k=0
Since hTk (x), Tm (x)i = αk δk,m , then when k = m, hTk (x), Tm (x)i = αk . So that hGi , Gj i = ˆ k,i = √αk Gk,i , then If we define G hGi , Gj i =
N +4 X
N +4 X
(3.43)
αk Gk,i Gk,j .
k=0
ˆ k,i G ˆ k,j = G
k=0
N +4 X
ˆ T )i,k G ˆ k,j . (G
(3.44)
k=0
In the same way, 2
2
hD Gi , D Gj i = =
*N +4 X
(M
2
G)T i,k Tk (x),
(M
2
G)T j,m Tm (x)
m=0
k=0
+4 N +4 N X X
N +4 X
+
(3.45)
2 T (M 2 G)T i,k (M G)j,m hTk (x), Tm (x)i
(3.46)
2 T (M 2 G)T i,k (M G)j,m αk δk,m .
(3.47)
k=0 m=0
=
+4 N +4 N X X
k=0 m=0
18
ˆ 2 = √αk M 2 , then and if we define M k,i k,i 2
N +4 X
2
hD Gi , D Gj i =
ˆ 2 G)T (M ˆ 2 G)T = (M i,k j,k
N +4 X
ˆ 2 G)T (M ˆ 2 G)k,j . (M i,k
(3.48)
k=0
k=0
ˆ 4 = √αk M 4 , we get that Now defining M k,i k,i *N +4 + N +4 N +4 X X X ˆ 4 G)T G ˆ hD 4 Gi , Gj i = (M (G)T = (M 4 G)T i,k k,j . j,m Tm (x) i,k Tk (x), m=0
k=0
(3.49)
k=0
Lastly, the inner product hk(x)Gi , Gj i will be calculated, k(x)Gi (x) = k(x)GT TN +4 =
N +4 X
k(x)Gk,i Tk (x).
(3.50)
k=0
If we substitute k(x) =
PN +4
m=0
km Tm (x) into the above equation then we get that
k(x)Gi (x) =
+4 N +4 N X X
km Gk,i Tm (x)Tk (x).
(3.51)
k=0 m=0
Since Tm (x)Tk (x) = 21 Tm+k (x) + 12 T|m−k| (x), then we may write the above equation as T
T
k(x)Gi (x) = [Wk gi ] TN +4 + [Wk∗ gi ] TN∗ +4 =
N +4 X
[Wk gi ]m Tm (x) +
N +4 X
∗ [Wk∗ gi ]m Tm+N +5 (x)
(3.52)
m=0
m=0
PN +4 where [Wk gi ]m = n=0 (Wk )m,n (gi )n . We will neglect the Wk∗ quantities since these terms will contribute nothing when taking the inner product with Gj (x). Since (gi )n = Gn,i , then k(x)Gi (x) =
+4 N +4 N X X
(Wk )m,n Gn,i Tm (x)
(3.53)
m=0 n=0
Upon taking the inner product of the above equation with Gj (x) we get that + *N +4 N +4 N +4 X X X Gp,j Tp (x) hk(x)Gi , Gj i = (Wk )m,n Gn,i Tm (x), p=0
m=0 n=0
=
+4 +4 N N +4 N X X X m=0 n=0 p=0
=
+4 N +4 N X X
(Wk )m,n Gn,i Gp,j hTm (x), Tp (x)i
αm (Wk )m,n Gn,i Gm,j .
(3.54)
m=0 n=0
ˆ k )m,n = √αm (Wk )m,n with G ˆ m,j = √αm Gm,j , the above equation may And thus after defining (W be expressed as hk(x)Gi , Gj i =
N +4 X
ˆ k G)T G ˆ (W i,m m,j .
(3.55)
m=0
Thus the left-hand side of equations (3.35a)-(3.35d) may be written in block matrix form as ˙ a(t) M11 0 M13 0 ˙ 0 M22 0 0 c(t) (3.56) M z˙ = ˙ , b(t) M31 0 M33 0 ˙ 0 0 0 M44 d(t) 19
where the matrix blocks are given by hD 2 Gi , D 2 Gj i hk(x)Gi , Gj i + , ri rj ri rj hk(x)Gi , Gj i , = ri rj
(M11 )j,i = p1 (M13 )j,i
(M22 )j,i = hGi , Gj i,
hk(x)Gi , Gj i = (M13 )j,i ri rj hD 2 Gi , D 2 Gj i hk(x)Gi , Gj i + , = p2 ri rj ri rj
(M31 )j,i = (M33 )j,i
(M44 )j,i = hGi , Gj i = (M22 )j,i . and the right-hand side of equations (3.35a)-(3.35d) 0 A12 0 A21 A22 A23 Az = 0 A32 0 A41 0 A43 where the matrix blocks are given by
may be written in block matrix form as a(t) A14 0 c(t) , (3.57) b(t) A34 d(t) A44
hD 2 Gi , D 2 Gj i hGi , Gj i + , rj rj hGi , Gj i =− , rj hD 4 Gi , Gj i hk(x)Gi , Gj i − , = −p1 ri ri = −δhk(x)Gi , Gj i,
(A12 )j,i = p1 (A14 )j,i (A21 )j,i (A22 )j,i
hk(x)Gi , Gj i , ri hGi , Gj i = (A14 )j,i , =− rj
(A23 )j,i = (A32 )j,i
hD 2 Gi , D 2 Gj i hGi , Gj i + , rj rj hk(x)Gi , Gj i = = (A23 )j,i , ri hD 4 Gi , Gj i hk(x)Gi , Gj i − , = −p2 ri ri = −δhk(x)Gi , Gj i = (A22 )j,i .
(A34 )j,i = p2 (A41 )j,i (A43 )j,i (A44 )j,i
Thus upon equating equations (3.56) and (3.57) we obtain the ential equations ˙ a(t) 0 M11 0 M13 0 0 ˙ c(t) A M 0 0 22 21 M z˙ = Az ⇒ ˙ = 0 M31 0 M33 0 b(t) ˙ A41 0 0 0 M44 d(t) with corresponding algebraic eigenvalue problem
ˆ λzˆ = M −1 Az. 20
following system of ordinary differA12 A22 A32 0
0 A23 0 A43
a(t) A14 0 c(t) A34 b(t) A44 d(t)
(3.58)
(3.59)
4
Numerical results
All the numerical experiments for equations (3.18) and (3.59) were performed with Matlab 2012b(64bit) on a Windows 7 system. According to the Matlab documentation, the command eig(D, Q) provides higher accuracy than eig(Q−1 D), and eig(Q−1 D) provides higher accuracy than does eig(Q\D).
4.1
Decoupled problem
When k(x) = 0, then the two beams equations (3.1a)-(3.1b) will be simplified to be two beam equations with exactly the same boundary conditions. In this part, we will consider the simplest case where p1 = p2 = 1, k(x) = 0, then (3.3a)-(3.3b) become ∂u1 = v1 ∂t ∂v1 ∂ 4 u1 =− ∂t ∂x4
(4.1a) (4.1b)
with boundary conditions (3.2a)-(3.2b), and ∂u2 = v2 ∂t ∂v2 ∂ 4 u2 =− ∂t ∂x4
(4.2a) (4.2b)
with boundary conditions (3.2c)-(3.2d). The system for u1 and u2 are thus decoupled into two beam equations with identical boundary conditions. In the following results from simulations, we calculate eigenvalues for N = 20 to N = 200 and only focus on those eigenvalues with magnitude less than 300. We present, in tabular form, part of the simulation results where N is a multiple of 20, starting with N = 20 up to N = 200. For a given N , we define that eigenvalue having the smallest magnitude as the first non-zero eigenvalue and that eigenvalue having magnitude closest to 300 as the last non-zero eigenvalue. Tables (4.1)-(4.3) are for the decoupled problem and we list the first and last non-zero eigenvalue as given by the Tau method [2] and when using L2 and Hilbert Space inner products to spatially discretize our system of partial differential equations. When using the Tau method, the variables u1 , u2 , v1 , v2 are expanded in the VN space [3] while for the L2 and Hilbert Space projection methods the variables are expanded in the VN +4 space. We note that with the L2 and Hilbert Space projection methods, the basis functions satisfy the boundary condition exactly while with the Tau method boundary conditions are only satisfied approximately. From Tables (4.1)-(4.3) we see that, regardless of whether we focus on the first non-zero eigenvalue or the last non-zero eigenvalue, the number of significant digits in the imaginary part of the eigenvalues given by the Tau method is only 5. But from results where the L2 and Hilbert Space inner products were used to spatially discretize our system of partial differential equation, we have at least 8 significant digits. The real part of the eigenvalues should be identically zero. Numerical results when using the L2 inner product as a projection method returned this exact value while with the Hilbert Space inner product method real parts were on the order of zero, to machine accuracy. The Tau method returned the poorest results for the real part of the eigenvalues.
4.2
Coupled problem 1 2
1 2
In this section, we choose k(x) = 0.1(e−16(x− 2 ) + e−16(x+ 2 ) ), p1 = p2 = 1, with the damping parameter δ = 0 or δ = 1. Tables (4.4)-(4.9) are for the coupled problem and we list the first and last non-zero eigenvalue as given by the Tau method and when using L2 and Hilbert Space inner
21
Table 4.1: Decoupled result of Tau method first non-zero eigenvalue last non-zero eigenvalue -6.705598014098E-14+2.467401094780 i 8.847124421951E-13+2.986688399691E+02 i 1.568884491585E-11+2.467402097159 i 1.761809974169E-11+2.985555331487E+02 i 7.165028143492E-11+2.467423720752 i -2.349985916899E-10+2.985555322398E+02 i -1.503783152077E-10+2.467613460975 i 7.668656776916E-10+2.985555313579E+02 i -2.032218010365E-10+2.466537807889 i -8.733368354941E-10+2.985555333780E+02 i 4.437809932871E-10+2.462436402624 i -3.840691884082E-09+2.985555342675E+02 i 3.377024176437E-09+2.476829981371 i 5.128603730863E-09+2.985556439165E+02 i 8.018927187781E-10+2.514746929390 i -7.461088354724E-09+2.985539417344E+02 i -1.022905982853E-09+2.131379055191 i -1.732687196718E-08+2.985547721941E+02 i -6.756858552262E-09+2.337543707013 i -2.889222089258E-08+2.985567659389E+02 i
N 20 40 60 80 100 120 140 160 180 200
Table 4.2: Decoupled result of L2 inner product method based on Chebyshev basis N first non-zero eigenvalue last non-zero eigenvalue 20 0+2.467401100269 i 0+2.986688399692E+02 i 40 0+2.467401100246 i 0+2.985555331330E+02 i 60 0+2.467401100172 i 0+2.985555331330E+02 i 80 0+2.467401100228 i 0+2.985555331330E+02 i 100 0+2.467401100475 i 0+2.985555331330E+02 i 120 0+2.467401100171 i 0+2.985555331330E+02 i 140 0+2.467401099878 i 0+2.985555331330E+02 i 160 0+2.467401099521 i 0+2.985555331330E+02 i 180 0+2.467401100908 i 0+2.985555320838E+02 i 200 0+2.467401100218 i 0+2.985555347055E+02 i
Table N 20 40 60 80 100 120 140 160 180 200
4.3: Decoupled result of Hilbert Space inner product method based on Chebyshev basis first non-zero eigenvalue last non-zero eigenvalue -1.306999380261E-14+2.467401100174 i -9.98326702631E-11+2.991678018448E+02 i -1.343450278756E-13+2.467401101420 i 1.169260017025E-11+2.985555331311E+02 i -7.554568084249E-14+2.467401072860 i 1.196110695876E-10+2.985555331703E+02 i 7.987145869450E-14+2.467401172821 i -1.103474043019E-10+2.985555330376E+02 i 2.083246001098E-15+2.467401191054 i 1.925940433517E-11+2.985555330213E+02 i -1.620447780703E-14+2.467401361637 i -8.172333694982E-11+2.985555328481E+02 i 1.418744016787E-12+2.467398981886 i -5.641281441202E-11+2.985555356969E+02 i 6.989363144672E-11+2.467403232593 i -1.007735885850E-10+2.985555305601E+02 i 9.026536089681E-10+2.467398107616 i -6.140341590391E-09+2.985555367307E+02 i 6.478229430914E-11+2.467397967167 i 7.649350863156E-10+2.985555373950E+02 i
22
products to spatially discretize our system of partial differential equations with δ = 0 and δ = 1. From Tables (4.4)-(4.9) we see that, regardless of whether we focus on the first non-zero eigenvalue or the last non-zero eigenvalue, the number of significant digits in the imaginary part of the eigenvalues given by the Tau method is only 5. But from results where the L2 and Hilbert Space inner products were used to spatially discretize our system of partial differential equation, we have at least 7 significant digits. The real part of the eigenvalues should be identically zero when δ = 0 and should be strictly negative when δ = 1. When δ = 0 the Hilbert Space projection method returns eigenvalues with real part identically zero and the L2 projection method returns eigenvalues with real part on the order of zero, to machine accuracy. The Tau method returns the poorest approximation to the real part of the eigenvalues and exhibits some stability problems relative to convergence of the imaginary part of the eigenvalues. When δ = 1 both the L2 and Hilbert Space projection method return eigenvalues with real part that have the desired property of being strictly negative and it would appear that the L2 projection method returns eigenvalues that converge more rapidly. This is not to be expected, given the superior behavior of the Hilbert space projection method in predicting eigenvalues when δ = 0 and warrants further investigation. The Tau method clearly returns the poorest approximation to the real part of the eigenvalues and also exhibits the same stability problems, relative to convergence of the imaginary part of the eigenvalues, that were present with δ = 0. Figures (4.1)-(4.4) show convergence behavior of the imaginary and real part of the last non-zero eigenvalue for the coupled problem. Figures (4.1) and (4.2) clearly exhibit convergence and stability problems that exist when using the Tau method. Figure (4.3) shows convergence behavior with the L2 projection method with δ = 1. This projection method appeared best when δ = 1. As mentioned above, this result is not expected, and warrants further investigation. Figure (4.4) shows
N 20 40 60 80 100 120 140 160 180 200
Table 4.4: Coupled result of Tau method, δ = 0 first non-zero eigenvalue last non-zero eigenvalue 1.813141231220E-13+2.467401113167 i -2.328331971668E-12+2.986688399691E+02 i 3.381363039473E-11+2.467399221474 i -2.609234060909E-11+2.985557343955E+02 i -6.585467406376E-10+2.467488951408 i 3.122536384275E-10+2.985557371651E+02 i 5.320960028680E-10+2.466886149922 i 7.871027188168E-10+2.985557410790E+02 i 1.433304326101E-10+2.467215116948 i 1.620106678028E-09+2.985557139647E+02 i -3.238582972387E-09+2.463253463769 i 4.838114339549E-09+2.985555452318E+02 i -1.207869623498E-08+2.453631441429 i -9.425007160857E-05+2.985558906559E+02 i 1.732980771025E-08+2.483283422922 i -2.332974594969E-08+2.985559961173E+02 i 7.815160878427E-02+2.277009910907 i -1.349500070194E-03+2.985560165261E+02 i 4.517035243871E-08+2.795036974019 i -7.191520240927E-08+2.985518060386E+02 i
N 20 40 60 80 100 120 140 160 180 200
Table 4.5: Coupled result of Tau method, δ = 1 first non-zero eigenvalue last non-zero eigenvalue -3.343092260546E-02+2.467268914505 i -4.012864610018E-02+2.986687939775E+02 i -3.046995364681E-02+2.467289557063 i -3.004363359992E-02+2.985557013780E+02 i -3.046591367662E-02+2.467276642906 i -3.046422238139E-02+2.985557041466E+02 i -3.046582954389E-02+2.467768941256 i -3.046440884390E-02+2.985557039390E+02 i -3.046453325873E-02+2.468520690905 i -3.046441030609E-02+2.985556991528E+02 i -3.047390197569E-02+2.470596691565 i -3.046442146404E-02+2.985557360775E+02 i -3.890404005919E-02+2.479766768709 i -3.037677448922E-02+2.985556384643E+02 i -3.055148980721E-02+2.434383459495 i -3.046503379262E-02+2.985553356671E+02 i -2.993809503398E-02+2.270872395463 i -3.046398072003E-02+2.985569115860E+02 i 8.985638519438E-02+2.379663589928 i -3.046351060327E-02+2.985539901649E+02 i
23
Table 4.6: Coupled result of L2 inner product method based on Chebyshev basis, δ = 0 N first non-zero eigenvalue last non-zero eigenvalue 20 4.794515479078E-14+2.467401100272 i 1.264266469292E-12+2.986688399691E+02 i 40 -2.345606348042E-12+2.467401100231 i -2.268518706217E-12+2.985555697252E+02 i 60 -3.055189299250E-11+2.467401099385 i -3.028798913063E-11+2.985555702384E+02 i 80 -1.640015196130E-11+2.467401096683 i -7.921982514425E-12+2.985555702456E+02 i 100 -6.000598263179E-11+2.467401112028 i 1.471674683159E-10+2.985555702282E+02 i 120 -1.439859495773E-10+2.467401119248 i -5.120891970016E-11+2.985555702114E+02 i 140 -1.185869050967E-10+2.467401069820 i -1.209905468348E-09+2.985555702877E+02 i 160 3.639652335283E-10+2.467401247723 i -5.078418258836E-10+2.985555706296E+02 i 180 3.667867737507E-11+2.467401227815 i -9.979618652723E-10+2.985555703877E+02 i 200 3.710481723475E-10+2.467400782526 i -4.916232485453E-09+2.985555692702E+02 i
Table 4.7: Coupled result of L2 inner product method based on Chebyshev basis, δ = 1 N first non-zero eigenvalue last non-zero eigenvalue 20 -6.077228184618E-03+2.472315451479 i -7.296147769658E-03+2.986688879395E+02 i 40 -5.539132572058E-03+2.471881225063 i -5.462420024517E-03+2.985555691914E+02 i 60 -5.538374080347E-03+2.471880612819 i -5.538890964086E-03+2.985555697036E+02 i 80 -5.538374021944E-03+2.471880614097 i -5.538924742400E-03+2.985555696935E+02 i 100 -5.538373921805E-03+2.471880622592 i -5.538924666957E-03+2.985555697245E+02 i 120 -5.538374235518E-03+2.471880655294 i -5.538925170105E-03+2.985555696625E+02 i 140 -5.538373989821E-03+2.471880557778 i -5.538924641325E-03+2.985555697754E+02 i 160 -5.538373246956E-03+2.471880631905 i -5.538924178500E-03+2.985555697666E+02 i 180 -5.538374384906E-03+2.471880740717 i -5.538927030167E-03+2.985555691089E+02 i 200 -5.538375905247E-03+2.471880313327 i -5.538926914480E-03+2.985555693472E+02 i
Table 4.8: Coupled result of Hilbert Space inner product method based on Chebyshev basis, δ = 0 N first non-zero eigenvalue last non-zero eigenvalue 20 0+2.467401100186 i 0+2.991653577748E+02 i 40 0+2.467401101144 i 0+2.985555331317E+02 i 60 0+2.467401120569 i 0+2.985555331102E+02 i 80 0+2.467401172647 i 0+2.985555329554E+02 i 100 0+2.467401436480 i 0+2.985555326753E+02 i 120 0+2.467399489321 i 0+2.985555351681E+02 i 140 0+2.467403285242 i 0+2.985555305551E+02 i 160 0+2.467409597270 i 0+2.985555225387E+02 i 180 0+2.467395532838 i 0+2.985555401061E+02 i 200 0+2.467413943868 i 0+2.985555050470E+02 i
24
Tau, TwoBeams,δ = 0 −4.7211297371507754e−12+298.6691086883846i 298.68 298.66
imaginary
298.64 298.62 298.6 298.58 298.56 298.54 20
40
60
80
100
40
60
80
100
120
140
160
180
200
120
140
160
180
200
−3
1
x 10
0.5
real
0
−0.5
−1
−1.5 20
terms
Figure 4.1: Last non-zero eigenvalue, coupled problem - Tau method, δ = 0 convergence behavior with the Hilbert Space projection method with δ = 0. This projection method clearly performed best with δ = 0.
4.3
Conclusion
From the comparison of the results for the decoupled and coupled problems, we can see that the L2 inner product and Hilbert Space inner product methods are more reliable than the Tau method since they provide better accuracy and better numerical stability. This does not mean the Tau method will perform poorly for every beam equation but only for the coupled beam equations presented in this paper. The L2 inner product method gives more promising results when the damping, δ = 1 while the Hilbert Space inner product method gives better results when the damping, δ = 0. It is possible that results from the Hilbert Space inner product method would perform better when δ = 1 if we
Table 4.9: Coupled result of Hilbert Space inner product method based on Chebyshev basis, δ = 1 N first non-zero eigenvalue last non-zero eigenvalue 20 -1.398357578269E-02+2.467363269590 i -8.813979442928E-03+2.991653559500E+02 i 40 -3.113953936625E-03+2.467399282939 i -5.164412664935E-03+2.985555327638E+02 i 60 -5.650443423560E-03+2.467394987625 i -5.456026119469E-03+2.985555328427E+02 i 80 -5.539505250050E-03+2.467395494831 i -5.546996595502E-03+2.985555342605E+02 i 100 -5.538671831835E-03+2.467395605785 i -5.558298756751E-03+2.985555337636E+02 i 120 -5.538901224057E-03+2.467393663824 i -5.533184301117E-03+2.985555393982E+02 i 140 -5.537829070998E-03+2.467396294386 i -5.499261838016E-03+2.985555265588E+02 i 160 -5.536659541852E-03+2.467400659306 i -5.603798305067E-03+2.985555462832E+02 i 180 -5.530276495927E-03+2.467407454090 i -5.271875342596E-03+2.985555327163E+02 i 200 -5.529889505552E-03+2.467407865463 i -5.662180738027E-03+2.985555355425E+02 i
25
Tau, TwoBeams,δ = 1 −0.040129225959461681+298.66906269470729i 298.68 298.66
imaginary
298.64 298.62 298.6 298.58 298.56 298.54 20
40
60
80
100
40
60
80
100
120
140
160
180
200
120
140
160
180
200
−0.024 −0.026 −0.028
real
−0.03 −0.032 −0.034 −0.036 −0.038 −0.04 −0.042 20
terms
Figure 4.2: Last non-zero eigenvalue, coupled problem - Tau method, δ = 1
/UseBC, TwoBeams, δ = 1 −0.0072961477696580901+298.66888793949005i 298.68 298.66
imaginary
298.64 298.62 298.6 298.58 298.56 298.54 20
40
60
80
100
40
60
80
100
120
140
160
180
200
120
140
160
180
200
−3
−4.5
x 10
−5
real
−5.5 −6 −6.5 −7 −7.5 20
terms
Figure 4.3: Last non-zero eigenvalue, coupled problem - L2 inner product method, δ = 1
26
H6UseBC, TwoBeams, δ = 0 0+299.20938085200459i 299.3 299.2 299.1
imaginary
299 298.9 298.8 298.7 298.6 298.5 298.4 20
40
60
80
100
40
60
80
100
120
140
160
180
200
120
140
160
180
200
−6
1
x 10
0 −1
real
−2 −3 −4 −5 −6 −7 −8 20
terms
Figure 4.4: Last non-zero eigenvalue, coupled problem - Hilbert Space inner product method, δ = 0 were able to more accurately approximate the scaling parameter, ri (see equation (A.6)). In this paper, we utilized the Tau method and both L2 and Hilbert Space inner product projection methods to reduce a system of partial differential equations to a system of ordinary differential equations. The algebraic eigenvalue problem associated with the system of ordinary differential equations was then solved using standard software. The eigenvalues for the system of ordinary differential equations are approximations to eigenvalues of the original system of partial differential equations and both the L2 and Hilbert Space inner product projection methods proved to be reasonable successful in approximating these eigenvalues. It will be valuable and instructive to investigate whether these methods can be applied to higher dimensional problems.
A
Appendix
Here are some properties of Chebyshev polynomials: T0 (x) = 1 T1 (x) = x Tn (x) = 2xTn−1 (x) − Tn−2 (x), n > 2
Tn (±1) = (±1)n , n = 0, 1, . . . Qk−1 2 2 dk (n+k) m=0 (n − m ) T (±1) = (±1) , k ≥ 1, n = 0, 1, . . . n k dx (2k − 1)!!
(A.1)
Particularly, d Tn (±1) = (±1)n+1 n2 dx 2 2 d2 n+2 n (n − 1) T (±1) = (±1) n dx2 3 27
(A.2)
Define inner product Z
hu, vi = then
1 −1
√
uv dx, 1 − x2
π, hTi (x), Tj (x)i = π2 , 0,
Suppose f (x) =
N X
an Tn (x),
n=0
i=j=0 i = j 6= 0 i 6= j
(A.3)
N X dk a(k) f (x) = n Tn (x), dxk n=0
(k)
the relationship between an and an is known as recurrence rule: (k)
(k)
cn−1 an−1 = 2na(k−1) + an+1 , n ≥ 1, k ≥ 1 n
(A.4)
where cj =
(
2, 1
j=0 j = 1, 2, . . .
PN (1) If we suppose aN = i=0 sn,i ai , n = 0, . . . , N then by substituting this sum into A.4 and equating coefficients of ai , it is fairly straightforward to show that the sn,i have the following properties: sn,i sn−1,n sn−1,i
= 0 , 2n = , cn−1 sn+1,i , = cn−1
n = 0, . . . , N, i = 0, . . . , n n = N, . . . , 1 n = N − 1, . . . , 1, i = n + 2, . . . , N .
Thus if we let M = [sn,i ] then from A.4 and mathematical induction, it can be shown that (k) a0 a0 (k) a1 a1 . = Mk .. . . . . (k) aN a
(A.5)
N
In order to normalize the maximum coefficient of u1 ,v1 ,u2 ,v2 in 3.19 to be approximately 1, we introduce the scale parameter γn through the following normalization, Z 1 2 (G′′n (x)) 1 √ dx = 1 γn2 −1 1 − x2 where Gn (x) is defined in equation (3.7). It is not possible to get a closed-form expression for γn . Therefore, we obtained approximate values to γn by performing a discrete least squares fit to the data generated from ri = a(1 + i)b = γi , i = 0, . . . , 25 . In this way we find that a ≈ 4.98479 and b ≈ 2.00304 and thus that ri = 4.98479(1 + i)2.00304 . which is still a good approximation to γi for i ≈ 200. 28
(A.6)
B B.1
Matlab Code main Listing B.1: main.m
1
% main f i l e
2 3 4 5
clear ; clear a l l ; close a l l ;
6 7 8 9
% N i s b e t w e e n n1 and n2 n1 = 2 0 ; n2 = 2 0 0 ;
10 11 12 13
% keep only e i g e n v a l u e s of l e n g t h < radius with % non−n e g a t i v e i m a g i n a r y p a r t radius = 300;
14 15 16 17 18
19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37
% %% %% Tau method : Decoupled problem k ( x ) =0, 1 beam f o r p = 1 %boundary c o n d i t i o n : 1 or 2 b u t i n t h i s problem 1 & 2 a r e same conditions f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = Tau1Beam ( n1 , r , p ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = Tau1Beam (m, r , p ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave1 ( 0 , r e s ( : ,m) , r , p , 0 , n1 , n2 ) ; end end end
38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57
%% Tau method : Coupled problem k ( x ) n o t e q u a l s t o 0 so t h a t 2 beams f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = Tau2Beam ( n1 , r ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = Tau2Beam (m, r ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave2 ( 0 , r e s ( : ,m) , r , 1 , n1 , n2 ) ; end end
58 59 60
%% L2 i n n e r p r o d u c t b a s e d on Chebshev b a s i s : %% v1 & v2 u s e boundary c o n d i t i o n : Decoupled problem (1 beam )
29
61
62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80
f o r p = 1 %boundary c o n d i t i o n : 1 or 2 ; b u t i n t h i s problem 1 & 2 a r e same conditions f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = L2UseBC1Beam ( n1 , r , p ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = L2UseBC1Beam (m, r , p ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave1 ( 1 , r e s ( : ,m) , r , p , 0 , n1 , n2 ) ; end end end
81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101
%% L2 i n n e r p r o d u c t b a s e d on Chebshev b a s i s : %% v1 & v2 u s e boundary c o n d i t i o n : c o u p l e d problem (2 beams ) f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = L2UseBC2Beam ( n1 , r ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = L2UseBC2Beam (m, r ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave2 ( 1 , r e s ( : ,m) , r , 1 , n1 , n2 ) ; end end
102 103 104 105 106
107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125
%% H e l b e r t Space I n n e r Product b a s e d on Chebshev b a s i s : Decoupled problem f o r p = 1 %boundary c o n d i t i o n : 1 or 2 ; b u t i n t h i s problem 1 & 2 a r e same conditions f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = HSUseBC1Beam ( n1 , r , p ) ;%% H Space tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = HSUseBC1Beam (m, r , p ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave1 ( 2 , r e s ( : ,m) , r , p , 0 , n1 , n2 ) ; end end end
126
30
127
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145
B.2
%% H e l b e r t Space I n n e r Product b a s e d on Chebshev b a s i s : Decoupled problem = 2 beams f o r r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = HSUseBC2Beam ( n1 , r ) ; tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f o r m = n1 +1: n2 tmp = HSUseBC2Beam (m, r ) ; pL = length ( r e s (m−n1 , : ) ) ; f o r k = 1 : pL r e s (m−n1 +1 ,k ) = G e t C l o s e s t P o i n t ( tmp , r e s (m−n1 , k ) ) ; end end f o r m = 1 : pL DrawSave2 ( 2 , r e s ( : ,m) , r , 1 , n1 , n2 ) ; end end
c Listing B.2: c.m
1 2 3 4 5 6 7 8
B.3
% c(k) function r = c ( k ) i f ( k==0) r = 2.0; else r = 1.0; end end
ChebyshevExpansion Listing B.3: ChebyshevExpansion.m
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18
B.4
% compute t h e c o e f f i c i e n t s o f t h e Chebyshev e x p a n s i o n s % f o r the given func function r = ChebyshevExpansion ( func , n ) r=zeros ( 1 , n+1) ; f1 = func (1) ; f 2 = f u n c ( −1) ; t =1:n−1; f o r m=0:n tmp = ( f 1 + f 2 ) / 2 ; tmp = tmp + f u n c ( cos ( pi ∗ t /n ) ) ∗ cos ( pi ∗m∗ t /n ) ’ ; i f (m==0 | | m==n ) r (m+1) = tmp/n ; else r (m+1) = 2∗tmp/n ; end f 2 = −f 2 ; end end
ChebyshevMultiplication Listing B.4: ChebyshevMultiplication.m
1 2 3 4 5 6
% f o r m u l a f o r f ( x ) g ( x ) where f ( x ) and g ( x ) a r e e x p a n s i o n s % i n term o f Chebyshev b a s i s function r = C h e b y s h e v M u l t i p l i c a t i o n ( c1 ) L = length ( c1 ) −1; r = zeros ( 2 ∗L+1 ,L+1) ; for n = 0 : L
31
for m = 0 : L r (m+n+1 , m+1) = r (m+n+1 , m+1) + c1 ( n+1) / 2 ; i f m>=n k = m − n; else k = n − m; end r ( k+1 , m+1) = r ( k+1 , m+1) + c1 ( n+1) / 2 ; end
7 8 9 10 11 12 13 14 15
end
16 17
end
B.5
GetClosestPoint Listing B.5: GetClosestPoint.m
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
% c h o o s e a complex number from ptmp which i s c l o s e s t t o a p o i n t function r = G e t C l o s e s t P o i n t ( ptmp , a p o i n t ) t = 0.01; l = 0; while l
