Use of Chebyshev Polynomials to Construct New Functions to Solve ...

Use of Chebyshev Polynomials to Construct New Functions to Solve Partial Differential Equations With High Accuracy

A Project Submitted to the Faculty of the Department of Mathematics and Statistics University of Minnesota, Duluth by

Qingtao Geng In partial fulfillment of the requirements for the degree of Master of Science in Applied and Computational Mathematics

Advisor: Steven A. Trogdon

September 9, 2010

Acknowledgements I am heartily thankful to my advisor, Dr. Steven A. Trogdon, whose encouragement, guidance, support, patience and carefulness from the initial to the final level enabled me to develop an understanding of the subject. Besides, I offer my regards and blessings to Dr. Zhuangyi Liu who helped me in any respect during the completion of the project. I would also like to thank my wife, Yu Guo, who typed in the initial draft of this paper and gave me much support during my study.

i

Abstract Based on Chebyshev polynomials, new basis functions, which satisfy the given boundary conditions, can be constructed to solve a system of PDE. With such basis functions, the boundary conditions are no longer a concern, and by exploiting the properties of Chebyshev polynomials, it is very easy to discretize the given PDE, and to solve the corresponding eigenvalue problem. It turns out that this technique gives very good accuracy. In this paper, the results given by this technique and Tau method are compared.

ii

Contents Acknowledgements

i

Abstract

ii

List of Tables

iv

List of Figures

v

List of Matlab Code

vi

1 Introduction

1

2 Notation and Formulae

3

3 Beam Equations 3.1 The Tau Method . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 New Chebyshev Basis . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Construction of the new basis . . . . . . . . . . . . . . . . . 3.2.2 Discretize with v1 and v2 satisfying boundary conditions . . 3.2.3 Discretize with v1 and v2 not satisfying boundary conditions

. . . . .

5 6 11 11 14 16

4 Numerical result 4.1 Separated problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Coupled problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

19 19 22 24

A Appendix

26

B Matlab Code B.1 report . . . . . . . . . . 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 NBUseBC1Beam . . . . B.11 NBUseBC2Beam . . . . B.12 NBNoBC1Beam . . . . . B.13 NBNoBC2Beam . . . . .

28 28 30 30 31 31 32 33 34 35 36 38 39 40

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

iii

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

List of Tables 4.1 4.2 4.3 4.4 4.5 4.6 4.7 4.8 4.9

Separated: result of (3.30) with (3.2a)-(3.2b) Separated: result of (3.36) with (3.2a)-(3.2b) Separated: result of (3.42) with (3.2a)-(3.2b) Separated: result of (3.30) with (3.2c)-(3.2d) Separated: result of (3.36) with (3.2c)-(3.2d) Separated: result of (3.42) with (3.2c)-(3.2d) Coupled: result of (3.30) . . . . . . . . . . . Coupled: result of (3.36) . . . . . . . . . . . Coupled: result of (3.42) . . . . . . . . . . .

iv

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

. . . . . . . . .

20 21 21 21 22 22 23 23 24

List of Figures 4.1 4.2 4.3

Eigenvalues given by (3.36) with N = 200 . . . . . . . . . . . . . . . Zoom-in of part of figure (4.1) . . . . . . . . . . . . . . . . . . . . . . Zoom-in of part of figure (4.1) . . . . . . . . . . . . . . . . . . . . . .

v

24 25 25

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

report.m . . . . . . . . . . . c.m . . . . . . . . . . . . . . ChebyshevExpansion.m . . . ChebyshevMultiplication.m GetClosestPoint.m . . . . . DrawSave1.m . . . . . . . . DrawSave2.m . . . . . . . . Tau1Beam.m . . . . . . . . Tau2Beam.m . . . . . . . . NBUseBC1Beam.m . . . . . NBUseBC2Beam.m . . . . . NBNoBC1Beam.m . . . . . NBNoBC2Beam.m . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

vi

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

. . . . . . . . . . . . .

28 30 30 31 31 32 33 34 35 36 38 39 40

1

Introduction

The spectral method is widely used to solve partial differential equations (PDEs). The method uses a linear combination of a set of orthogonal basis functions to approximate the solution to the PDE and employs a projection procedure to obtain a system of equations for the coefficients multiplying the basis functions. 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 √ hu, vi = dx . 1 − x2 −1 In general, for a PDE of order 2k Lu = g, −1 ≤ x ≤ 1

(1.1)

Bl,i u(−1) = 0, Br,i u(1) = 0, i = 1, . . . , k

(1.2)

with boundary conditions

where L, Bl,i and 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)

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.3)

which is equivalent to N X n=0

an hLTn (x), Ti (x)i = hg, Ti (x)i, i = 0, . . . , N .

(1.4)

To satisfy the boundary conditions (1.2), we require that Bl,i uN (−1) = 0, Br,i uN (1) = 0, i = 1, . . . , k

(1.5)

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 .

n=0

1

(1.6)

In equations (1.4) and (1.6), there are N + 1 + 2k equations with N + 1 unknowns. To make number of equations and unknowns consistent, the inner product in (1.4) is only applied for i = 0, . . . , N − 2k; the remaining 2k equations in (1.4) are dropped from the analysis. These dropped 2k equations may be used to define tau coefficients through τj =

N X n=0

an hLTn (x), TN −j (x)i − hg, TN −j (x)i, j = 0, . . . , 2k − 1 .

(1.7)

The tau coefficients are measures of how accurate is the approximation, uN . They may be calculated after solving equations (1.4) for i = 0, . . . , N − 2k with boundary conditions (1.6) for coefficients an . The closer the tau coefficients are to zero the greater the accuracy of the approximation, uN . This procedure of dropping equations from equations like (1.4) to accommodate boundary conditions like (1.6) and then using tau-coefficient equations like (1.7) to assess the accuracy of an approximation is referred to as the tau method. That which has been said above relative to PDEs may also be applied to eigenvalue problems; the only difference being that instead of determining approximating coefficients, an , approximations to eigenvalues and eigenvectors are determined. The tau coefficients are then used to assess the accuracy of these eigenvalues and eigenvectors. One may now ask the question, is there any way to impose the boundary conditions (1.6) on uN while making use of all the N + 1 equations in (1.4)? The answer is, yes there is. Chebyshev polynomials may be used to construct new basis functions which satisfy boundary conditions (1.6). In (1.6), there are a total of 2k boundary conditions. If we let Sn (x) = Tn (x) +

2k X j=1

sn,j Tn+j (x), sn,j ∈ R, n = 0, 1, . . .

(1.8)

then, upon substituting (1.8) into (1.6) we obtain Bl,i Tn (−1) +

2k X

sn,j Bl,i Tn+j (−1) = 0, i = 1, . . . , k

(1.9)

sn,j Br,i Tn+j (1) = 0, i = 1, . . . , k .

(1.10)

j=1

Br,i Tn (1) +

2k X j=1

There are 2k unknowns sn,j , j = 1, . . . , 2k, and 2k algebraic equations in (1.9) and (1.10). If this system of equations is solvable the values of sn,j can be determined and thus the new basis functions Sn (x) can be constructed. Since each Sn (x) satisfies all of the boundary conditions in (1.6), u can be approximated by uN =

N X n=0

2

an Sn (x),

and (1.4) can be solved directly as à ! 2k N X X sn,j hLTn+j (x), Ti (x)i = hg, Ti (x)i, i = 0, . . . , N . an hLTn (x), Ti (x)i + n=0

j=1

(1.11)

2

Notation and Formulae

For convenience, we introduce some notation. Suppose f (x) =

N X

fi Ti (x) = f T TN

i=0

where

   T0 (x) f0     .. f =  ...  , TN =   . 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 (x) −→ f

where VN is the Chebyshev function space defined in Section 1. This notation 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 (x) −→ f ⇒

dk VN f (x) −→ M kf , k ≥ 1 k dx

(2.1)

where the matrix, M is of order (N + 1) × (N + 1) and is defined in Appendix A. Now let   g0 VN   g(x) −→ g =  ...  . gN Then, using the multiplication law,

1 1 Ti (x)Tj (x) = Ti+j (x) + T|i−j| (x) 2 2

3

for Chebyshev polynomials the product f (x)g(x) may be approximated as, f (x)g(x) =

N X

fi Ti (x)

i=0

=

N X

gj Tj (x)

j=0

N X N X

fi gj Ti (x)Tj (x)

i=0 j=0

N N £ ¤ 1 XX = fi gj Ti+j (x) + T|i−j| (x) 2 i=0 j=0 2N

N

1X 1X S1,i Ti (x) + S2,i Ti (x) = 2 i=0 2 i=0

where

S1,i =

X

fp gq , 0 ≤ i ≤ 2N

(2.2)

X

fp gq , 0 ≤ i ≤ N .

(2.3)

p+q=i 0≤p,q≤N

and S2,i =

|p−q|=i 0≤p,q≤N

From (2.2) we have that S1,0 = f0 g0 S1,1 = f0 g1 + f1 g0 S1,2 = f0 g2 + f1 g1 + f2 g0 .. . S1,N = f0 gN + f1 gN −1 + . . . + fN g0 and thus that     

S1,0 S1,1 .. . S1,N





    =  

f0 f1 .. .



f0

fN fN −1 . . . f1 f0

Likewise, from (2.3) we have that

   

g0 .. . .. . gN



   ≡ W1 g . 

S2,0 = f0 g0 + f1 g1 + . . . + fN gN S2,1 = f0 g1 + f1 g2 + . . . + fN −1 gN +f1 g0 + f2 g1 + . . . + fN gN −1 .. . S2,N = f0 gN + fN g0 4

and thus that    S2,0  S2,1       ..  =   .   S2,N

f0 f1 .. . fN

Thus, we get that

 f1 f2 ... fN −1 fN  f0 + f1 f1 + f2 . . . fN −2 + fN fN −1     0 ... 0 f0 2N

g0 .. . .. . gN



   ≡ W2 g . 

N

1X 1X S1,i Ti (x) + S2,i Ti (x) f (x)g(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 [W1 g]T TN + [W2 g]T TN + [W3 g]T TN∗ 2 2 2 · ¸T ¸T · 1 1 = (W1 + W2 )g TN + W3 g TN∗ 2 2 =

where

 TN +1 (x)   .. TN∗ =   . . T2N (x) 

If we let Wf = 21 (W1 + W2 ), and Wf∗ = 12 W3 , then

£ ¤T f (x)g(x) = [Wf g]T TN + Wf∗ g TN∗

(2.4)

where Wf∗ is not given explicitly because it consists of only terms of order higher than N , which become zero after applying inner product with Tk (x), k = 0, 1, . . . , N .

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 = −p − k(x) (u1 − u2 ) − δk(x) 1 2 4 ∂t ∂x ∂t

(3.1a)

∂ 2 u2 ∂u2 ∂ 4 u2 = −p + k(x) (u − u ) − δk(x) 2 1 2 ∂t2 ∂x4 ∂t

(3.1b)

5

with boundary conditions ∂ 2 u1 ∂ 2 u1 (t, −1) = (t, 1) = 0 ∂x2 ∂x2

(3.2a)

∂ 3 u1 ∂ 3 u1 (t, −1) = (t, 1) = 0 ∂x3 ∂x3

(3.2b)

u2 (t, 1) = 0

(3.2c)

∂u2 ∂u2 (t, −1) = (t, 1) = 0 ∂x ∂x

(3.2d)

u2 (t, −1) =

where p1 > 0, p2 > 0, and δ ≥ 0 are given parameters and k(x) is a known function. 1 2 If we introduce the substitutions, v1 = ∂u and v2 = ∂u , the original equations ∂t ∂t (3.1a) and (3.1b) become ∂u1 ∂t ∂v1 ∂t ∂u2 ∂t ∂v2 ∂t

=

(3.3a)

v1

∂ 4 u1 = −p1 4 − k(x)(u1 − u2 ) − δk(x)v1 ∂x

(3.3b)

=

(3.3c)

v2

= −p2

∂ 4 u2 + k(x)(u1 − u2 ) − δk(x)v2 . ∂x4

(3.3d)

Section 3.1 describes how to apply the Tau method to the above equations, section 3.2 describes how to construct new basis functions to solve the above equations, and section 4 compares the numerical results given by these two techniques.

3.1

The Tau Method

To develop the discretized system of equations associated with the tau method, we first expand the variables appearing in equations (3.3a) - (3.3d) in the function space VN as,   ( a0 (t) ∂u1 VN ˙ −→ a(t) VN  ..  u1 −→ a(t) =  .  ⇒ ∂∂t 4u V N 1 −→ M 4 a(t) ∂x4 aN (t)   ( b0 (t) ∂u2 VN ˙ −→ b(t) VN   u2 −→ b(t) =  ...  ⇒ ∂∂t 4u VN 4 2 4 −→ M b(t) ∂x bN (t)   c0 (t) ∂v1 VN VN   ˙ −→ c(t) v1 −→ c(t) =  ...  ⇒ ∂t cN (t) 6

 d0 (t) ∂v2 VN ˙ VN   v2 −→ d(t) =  ...  ⇒ −→ d(t) ∂t dN (t) 

and

 k0 VN   k(x) −→ k =  ...  , kN 

where the superposed dot denotes time differentiation. Next, for an arbitrary vector r(t) in VN , we introduce the partitioning:     rN −3 (t) r0 (t) ¸ · rf (t)     .. .. , rf (t) =  r(t) =   , rl (t) =  . . rl (t) rN (t) rN −4 (t)

and define the vectors R1 and R−1 through    T0 (1) T0 (−1)  T1 (1)   T1 (−1)    R1 =   , R−1 =  .. ..    . . TN (1) TN (−1)

(3.4)



  . 

The above definitions, partitioning and function space expansions are next applied to the boundary conditions, (3.2a) - (3.2d). Using the expansion for u1 and the vectors R1 and R−1 , we write the boundary conditions, (3.2a)-(3.2b), as £ 2 ¤T M a(t) £ 2 ¤T M a(t) £ 3 ¤T M a(t) £ 3 ¤T M a(t)

R−1 R1 R−1 R1

£ ¤T = 0 ⇒ [a(t)]T M 2 £ ¤T = 0 ⇒ [a(t)]T M 2 £ ¤T = 0 ⇒ [a(t)]T M 3 £ ¤T = 0 ⇒ [a(t)]T M 3 T

= 0 ≡ [a(t)]T S1 = 0

R−1

= 0 ≡ [a(t)]T S2 = 0

R1

= 0 ≡ [a(t)]T S3 = 0

R−1

= 0 ≡ [a(t)]T S4 = 0 , (3.5)

R1 T

where S1 and S3 are respectively, [M 2 ] and [M 3 ] applied to R−1 and S2 and S4 T T are respectively, [M 2 ] and [M 3 ] applied to R1 . Using the partitioning of (3.4), we can write equation (3.5) in partitioned form as, ¸ h i· S T T 1,f S2,f S3,f S4,f =0 [af (t)] , [al (t)] S1,l S2,l S3,l S4,l or, alternatively as

[af (t)]T P1 + [al (t)]T P2 = 0 (3.6) ¤ ¤ £ £ where P1 = S1,f S2,f S3,f S4,f and P2 = S1,l S2,l S3,l S4,l . Assuming P2 is invertible, we can solve (3.6) for al (t), 7

£ ¤T al (t) = − P1 P2−1 af (t).

(3.7)

As in the development above that led to equations (3.6) and (3.7), we now use the expansion for u2 and write the boundary conditions, (3.2c)-(3.2d), as [b(t)]T R−1 [b(t)]T R1 [M b(t)]T R−1 [M b(t)]T R1

=0 ⇒

=0 ⇒

[b(t)]T M T R−1 [b(t)]T M T R1

= 0 ≡ [b(t)]T S5 = 0

= 0 ≡ [b(t)]T S6 = 0

= 0 ≡ [b(t)]T S7 = 0

= 0 ≡ [b(t)]T S8 = 0

(3.8)

which may be written in partitioned form as ¸ i· S h T T 5,f S6,f S7,f S8,f =0 [bf (t)] , [bl (t)] S5,l S6,l S7,l S8,l

or, alternatively as

[bf (t)]T P3 + [bl (t)]T P4 = 0 (3.9) ¤ ¤ £ £ where P3 = S5,f S6,f S7,f S8,f and P4 = S5,l S6,l S7,l S8,l . Assuming P4 is invertible, we can solve (3.9) for bl (t), £ ¤T bl (t) = − P3 P4−1 bf (t).

(3.10)

£ ¤T b˙ l (t) = − P3 P4−1 b˙ f (t).

(3.12)

We note that P1 and P3 are of order (N − 3) × 4 and that P2 and P4 are of order 4 × 4. From (3.7) and (3.10), we see that £ ¤T a˙ l (t) = − P1 P2−1 a˙ f (t) (3.11)

and

Furthermore, using equations (3.7) and (3.10), we may eliminate the last four elements from the vectors a(t) and b(t) and represent those vectors as ¸ · IN −3 ¤T af (t) ≡ C1 af (t) £ (3.13) a(t) = − P1 P2−1 ¸ · IN −3 ¤T bf (t) ≡ C2 bf (t), £ (3.14) b(t) = − P3 P4−1

where the partitioned matrices C1 and C2 are of order (N + 1) × (N − 3) and IN −3 denotes the (N − 3) × (N − 3) identity matrix. Thus far we have discretized and accounted for the boundary conditions, (3.2a) - (3.2d). The next step in obtaining the system of equations associated with the tau method is to discretize the system of differential equations, (3.3a) - (3.3d). Using the function space expansions for u1 and v1 in VN , we obtain from equation (3.3a) that, T ˙ [a(t)] TN = [c(t)]T TN .

8

Upon taking the inner product of the above with Ti (x), i = 0, . . . , N , we get that, a˙ i (t) = ci (t), i = 0, . . . , N, and thus, with the partitioning of (3.4), that a˙ f (t) = cf (t)

(3.15)

a˙ l (t) = cl (t).

(3.16)

and Equations (3.15) and (3.16) with equation (3.11) imply that £ ¤T cl (t) = − P1 P2−1 cf (t)

and hence, as in equation (3.13), that

(3.17)

c(t) = C1 cf (t).

Similarly, upon using the function space expansions for u2 and v2 in VN , we obtain from (3.3c), that b˙ f (t) = df (t) (3.18) and

b˙ l (t) = dl (t)

(3.19)

and thus, with equation (3.12), that £ ¤T dl (t) = − P3 P4−1 df (t).

Hence, as in equation (3.14), we get that

(3.20)

d(t) = C2 df (t).

Equation (3.17) means that v1 must satisfy the same boundary conditions, (3.2a)(3.2b), as u1 and equation (3.20) means that v2 must satisfy the same boundary conditions, (3.2c)-(3.2d), as u2 . In fact, if u1 and u2 are sufficiently smooth, then it can be shown from (3.3a) and (3.3c) that ∂ 2 v1 ∂ ∂ 2 u1 ∂ 2 ∂u1 (t, ±1) = (t, ±1) = (t, ±1) ∂x2 ∂x2 ∂t ∂t ∂x2 ∂ 3 ∂u1 ∂ ∂ 3 u1 ∂ 3 v1 (t, ±1) = (t, ±1) = (t, ±1) ∂x3 ∂x3 ∂t ∂t ∂x3 ∂u2 v2 (t, ±1) = (t, ±1) ∂t ∂v2 ∂ ∂u2 ∂ ∂u2 (t, ±1) = (t, ±1) = (t, ±1) ∂x ∂x ∂t ∂t ∂x

=0

(3.21)

=0

(3.22)

=0

(3.23)

= 0.

(3.24)

which are the same conditions as (3.2a) - (3.2d), but for v1 and v2 . 9

The remaining equations to be discretized are equations (3.3b) and (3.3d). Using (2.4) for functions expanded in the function space, VN we get that k(x)u1 k(x)u2 k(x)v1 k(x)v2

= = = =

[Wk a(t)]T TN + [Wk∗ a(t)]T TN∗ [Wk b(t)]T TN + [Wk∗ b(t)]T TN∗ [Wk c(t)]T TN + [Wk∗ c(t)]T TN∗ [Wk d(t)]T TN + [Wk∗ d(t)]T TN∗ .

With the function space expansion of u1 , u2 and v1 in VN , we then obtain from (3.3b) that, £ ¤T ˙ T TN = −p1 M 4 a(t) TN − [Wk a(t)]T TN − [Wk∗ a(t)]T TN∗ [c(t)] + [Wk b(t)]T TN + [Wk∗ b(t)]T TN∗ −δ [Wk c(t)]T TN − δ [Wk∗ c(t)]T TN∗

Upon taking the inner product of the above with Ti (x), i = 0, . . . , N − 4, we get that £ ¤ (3.25) c˙i (t) = −p1 M 4 i a(t) − [Wk]i a(t) + [Wk]i b(t) − δ [Wk]i c(t), where [ ∗ ]i denotes the ith row of the indicated matrix. Equations (3.25) may be written in matrix-vector form as,

c˙ f (t) = A1 a(t) − Ba(t) + Bb(t) − δBc(t)

(3.26)

where A1 , of order (N − 3) × (N + 1), is the submatrix obtained from the first N − 3 rows of −p1 M 4 and B, of order (N − 3) × (N + 1), is the submatrix obtained from the first N − 3 rows of Wk. Using equations (3.13), (3.14) and (3.17), we may simplify equation (3.26) as c˙ f (t) = (A1 − B)C1 af (t) + BC2 bf (t) − δBC1 cf (t),

(3.27)

where only the first N − 3 elements of vectors appear as unknowns. Similarly, by substituting the function space expansion for u1 , u2 and v2 into equation (3.3d), taking the inner product of the resulting equation with Ti (x), i = 0, . . . , N − 4 and using equations (3.13), (3.14) and (3.20), we may deduce that d˙f (t) = BC1 af (t) + (A2 − B)C2 bf (t) − δBC2 df (t),

(3.28)

where A2 , of order (N − 3) × (N + 1), is the submatrix obtained from the first N − 3 rows of −p2 M 4 . If we now let   af (t)  bf (t)   z(t) =   cf (t)  df (t) 10

then equations (3.15), (3.18), (3.27) and (3.28) may be combined to give   0 0 IN −3 0  0 0 0 IN −3   z(t)  ˙ = z(t)   (A1 − B)C1 BC2 −δBC1 0 BC1 (A2 − B)C2 0 −δBC2 ≡

(3.29)

Dz(t).

The above equations are the discretized system of equations associated with the tau method. The eigenvalue problem for these equations is obtained by setting z(t) = ˆ where eλt z,           ˆ a d0 c0 b0 a0 ˆ   b .  ˆ =  ..  , cˆ =  ..  , dˆ =  ..  .  ˆ= zˆ =   .   .   .   ..  , b  cˆ  , a dN −4 cN −4 bN −4 aN −4 dˆ

Therefore, since

λeλt zˆ = Deλt zˆ we obtain D zˆ = λzˆ

(3.30)

which is the eigenvalue problem corresponding to the tau equations, (3.29).

3.2

New Chebyshev Basis

To demonstrate the new basis technique presented in section 1, we construct the new function spaces GN and HN first, where each basis of GN satisfies boundary conditions (3.2a)-(3.2b) and each basis of HN satisfies boundary conditions (3.2c)-(3.2d). This is presented in section 3.2.1. Then we expand u1 and v1 in GN , u2 and v2 in HN because they should satisfy the associated boundary conditions, and solve (3.3a)-(3.3d) with boundary conditions (3.2a)-(3.2d). This is presented in section 3.2.2. However, we also tried expanding u1 in GN , u2 in HN , but v1 and v2 in VN , which means v1 and v2 do not satisfy the given boundary conditions. We used this expansions to solve the given equations, and compared the significant digits. This is presented in section 3.2.3. 3.2.1

Construction of the new basis

To approximate u1 , a subspace of VN +4 , in which all the basis functions satisfy boundary conditions (3.2a) and (3.2b), is constructed as below.

11

Let GN = span {Gn (x), n = 0, . . . , N } where Gn (x) = Tn (x) + a1 Tn+1 (x) + a2 Tn+2 (x) + a3 Tn+3 (x) + a4 Tn+4 (x) 4 X = Tn (x) + ai Tn+i (x) i=1

Plugging the above Gn (x) into (3.2a) and (3.2b), and using the formula (A.1) from Appendix A, we get 4

2 2 n2 (n2 − 1) X d2 i (n + i) [(n + i) − 1] ai G (−1) = 0 ⇒ + =0 (−1) n dx2 3 3 i=1 4

d2 n2 (n2 − 1) X (n + i)2 [(n + i)2 − 1] ai G (1) = 0 ⇒ + =0 n dx2 3 3 i=1 4

2 2 2 d3 n2 (n2 − 1)(n2 − 4) X i (n + i) [(n + i) − 1] [(n + i) − 4] ai G (−1) = 0 ⇒ + =0 (−1) n dx3 15 15 i=1 4

d3 n2 (n2 − 1)(n2 − 4) X (n + i)2 [(n + i)2 − 1] [(n + i)2 − 4] ai Gn (1) = 0 ⇒ + =0 dx3 15 15 i=1

These are algebraic equations of four unknowns, ai , i = 1, . . . , 4, we solve them with Mathematica and get a1 = a3 = 0, a2 (n) =

(n − 1)n2 (n + 1)2 −2(n − 1)n2 , a (n) = 4 (n + 2)(n + 3)2 (n + 5)(n2 + 7n + 12)2

therefore Gn (x) = Tn (x) + a2 (n)Tn+2 (x) + a4 (n)Tn+4 (x) It can be shown that G0 (x), G1 (x), . . . , GN (x) are linearly independent, and they span the subspace, GN , of VN +4 . Similarly a subspace HN = span {Hn (x), n = 0, . . . , N } of VN +4 , which satisfies boundary conditions (3.2c) P4and (3.2d), can be constructed as below. Let Hn (x) = Tn (x) + i=1 bi Tn+i (x), then plug Hn (x) into (3.2c) and (3.2d), and use (A.1), we can get Hn (−1) = 0 ⇒ 1 + Hn (1) = 0 ⇒ 1 +

4 X

i=1 4 X

bi = 0

i=1 4 X

d Hn (−1) = 0 ⇒ n2 + dx d Hn (1) = 0 ⇒ n2 + dx 12

(−1)i bi = 0

(−1)i (n + i)2 bi = 0

i=1

4 X i=1

(n + i)2 bi = 0

Likewise we can solve these algebraic equations for the unknowns bi , i = 1, . . . , 4 with Mathematica, and get b1 = b3 = 0, b2 (n) =

n+1 −2(n + 2) , b4 (n) = n+3 n+3

therefore Hn (x) = Tn (x) + b2 (n)Tn+2 (x) + b4 (n)Tn+4 (x) It can be shown that H0 (x), H1 (x), . . . , HN (x) are linearly independent, and they span the subspace, HN , of VN +4 . Easy to see, every g(x) ∈ GN must be representable in VN +4 . In fact, we have VN +4

G

where

N g(x) −→ g ⇒ g(x) −→ Gg



1

0

0

...

 1 0  0   a (0) 0 1  2   0 a2 (1) 0   a4 (0) 0 a2 (2) G=  0 a4 (1) 0   ..  . 0 a4 (2)  . ..  .. . 0   . . .. ..  .. . 0 0 0

0 .. . .. . .. .



        0   1    0   a2 (N )    0  a4 (N )

has order (N + 5) × (N + 1), and represents the linear transformation from GN to VN +4 . Similarly it can be shown that VN +4

H

where

N h(x) −→ h ⇒ h(x) −→ Hh



1

0

0

 1 0  0   b (0) 0 1  2   0 b2 (1) 0   b4 (0) 0 b2 (2) H=  0 b (1) 0 4   ..  . 0 b4 (2)  . ..  .. . 0   . .. . ..  .. . 0 0 0

...

0 .. . .. . .. .



        0   1    0   b2 (N )    0  b4 (N )

has order (N + 5) × (N + 1), and represents the linear transformation from HN to VN +4 . 13

3.2.2

Discretize with v1 and v2 satisfying boundary conditions

As shown in section 3.1, v1 should satisfy the boundary conditions of u1 , and v2 should satisfy the boundary conditions of u2 , therefore u1 and v1 are expanded in GN , and u2 and v2 are expanded in HN . Furthermore, u1 , u2 , v1 and v2 should be mapped back into VN +4 because we need to apply inner product to them with Tn (x), n = 0, . . . , N , thus k(x) is expanded in VN +4 since it will multiply with functions from VN +4 instead of VN . The expansions are as following.  V   N +4   a0 (t) u1 −→ Ga(t) GN   VN +4 1 u1 −→ a(t) =  ...  ⇒ ∂u ˙ −→ Ga(t) ∂t   V 4 N +4  ∂ u1 aN (t) −→ M 4 Ga(t) ∂x4  V   N +4   b0 (t) u2 −→ Hb(t) HN   VN +4 2 ˙ u2 −→ b(t) =  ...  ⇒ ∂u −→ H b(t) ∂t    ∂ 4 u2 VN +4 bN (t) −→ M 4 Hb(t) ∂x4   ( VN +4 c0 (t) v1 −→ Gc(t) GN  ..  v1 −→ c(t) =  .  ⇒ ∂v1 VN +4 ˙ −→ Gc(t) ∂t cN (t)   ( VN +4 d0 (t) v2 −→ Hd(t) HN   v2 −→ d(t) =  ...  ⇒ ∂v2 VN +4 ˙ −→ H d(t) ∂t dN (t)   k0 VN +4   k(x) −→ k =  ...  kN +4 From (3.3a), we get

N X

a˙i (t)Gi (x) =

i=0

N X

ci (t)Gi (x)

i=0

Since Gi (x), i = 0, . . . , N are linearly independent, so 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.31)

where IN +1 is identity matrix of order N + 1. This is the discretized equation associated with (3.3a). 14

Similarly, we can get the discretized equation associated with (3.3c) as ˙ IN +1 b(t) = IN +1 d(t)

(3.32)

The remaining equations to be discretized are equations (3.3b) and (3.3d). Using (2.4) for functions expanded in the function space, VN +4 we get that k(x)u1 = [WkGa(t)]T TN +4 + [Wk∗ Ga(t)]T TN∗ +4 k(x)u2 = [WkHb(t)]T TN +4 + [Wk∗ Hb(t)]T TN∗ +4 k(x)v1 = [WkGc(t)]T TN +4 + [Wk∗ Gc(t)]T TN∗ +4 k(x)v2 = [WkHd(t)]T TN +4 + [Wk∗ Hd(t)]T TN∗ +4 With the function space expansion of u1 , u2 and v1 in VN +4 , we then obtain from (3.3b) that, ˙ T TN +4 = [Gc(t)]

£

¤T −p1 M 4 Ga(t) TN +4 − [WkGa(t)]T TN +4

− [Wk∗ Ga(t)]T TN∗ +4 + [WkHb(t)]T TN +4

+ [Wk∗ Hb(t)]T TN∗ +4 − δ [WkGc(t)]T TN +4 −δ [Wk∗ Gc(t)]T TN∗ +4

There are totally N + 5 equations after taking inner product with Tj (x), j = 0, . . . , N + 4, here we choose the first N + 1 equations and drop the last 4. In fact, if k(x) is 0, i.e. the two beams are independent from each other, the right hand side of the last 4 equations become 0 which makes the system of equations inconsistent. Upon taking inner product with Tj (x), j = 0, . . . , N , we get £ ¤ ˙ = −p1 M 4 G j a(t) − [WkG]j a(t) + [WkH]j b(t) [G]j c(t) −δ [WkG]j c(t),

j = 0, . . . , N

These can be written in matrix-vector form as, ˙ = A1 a(t) − B1 a(t) + B2 b(t) − δB1 c(t) G1 c(t)

(3.33)

This is the discretized equation corresponding to (3.3b), where G1 , of order (N + 1) × (N + 1), is a submatrix obtained from the first N + 1 rows of G, A1 , of order (N + 1) × (N + 1) is a submatrix obtained from the first N + 1 rows of −p1 M 4 G, B1 , of order (N + 1) × (N + 1) is a submatrix obtained from the first N + 1 rows of WkG, and B2 , of order (N + 1) × (N + 1), is a submatrix obtained from the first N + 1 rows of WkH. Similarly substituting the function space expansion of u1 , u2 and v2 into (3.3d), taking inner product with Tj (x), j = 0, . . . , N , we can get £ ¤ ˙ [H]j d(t) = −p2 M 4 H j b(t) + [WkG]j a(t) − [WkH]j b(t) −δ [WkH]j d(t),

15

j = 0, . . . , N

These can be written in matrix-vector form as, ˙ = A2 b(t) + B1 a(t) − B2 b(t) − δB2 d(t) H1 d(t)

(3.34)

This is the discretized equation corresponding to (3.3d), where A2 , of order (N + 1) × (N + 1), is the submatrix obtained from the first N + 1 rows of −p2 M 4 H. If we now let   a(t)  b(t)   z(t) =   c(t)  d(t)

then (3.31)-(3.34) may be  IN +1 0 0  0 I 0 N +1   0 0 G1 0 0 0 ˙ = Dz(t) ≡ Qz(t)

combined to give   0 0 0 IN +1 0  0  0 0 0 I N +1  z(t) ˙ =  A 1 − B1 0  B2 −δB1 0 H1 B1 A 2 − B2 0 −δB2



  z(t)  (3.35)

The above equations are the discretized system of equations associated with the new basis technique while v1 and v2 satisfying their boundary conditions. The eigenˆ where value problem for these equations is obtained by setting z(t) = eλt z,           ˆ a d0 c0 b0 a0 ˆ   b .  ˆ =  ..  , cˆ =  ..  , dˆ =  ..   ˆ= zˆ =   .   .   .   ..  , b  cˆ  , a dN cN bN aN dˆ therefore we obtain

λQeλt zˆ = Deλt zˆ ⇒ Q−1 D zˆ = λzˆ

(3.36)

which is the eigenvalue problem corresponding to the equation (3.35). 3.2.3

Discretize with v1 and v2 not satisfying boundary conditions

As an experiment, the procedure in section 3.2.2 is repeated without applying boundary conditions to v1 and v2 , i.e. v1 and v2 are expanded in VN while u1 and u2 are expanded in GN and HN , respectively, therefore k(x) should be expanded as k1 in VN +4 , and k2 in VN because it will multiply with functions from VN +4 and VN as well. The twice expansions of k(x) is caused by (2.4), which defines the multiplication between two functions from the same function space, VN . If we develop another formula for the multiplication between one function from VN and another function from VN +4 , then k(x) can be expended only once in VN . Now we choose expansions as below.  V   N +4   a0 (t) u1 −→ Ga(t) GN   VN +4 1 a(t) =  ...  ⇒ ∂u u1 −→ ˙ −→ Ga(t) ∂t   4 u VN +4  ∂ aN (t) 1 −→ M 4 Ga(t) ∂x4

16

 V  N +4   b0 (t) u2 −→ Hb(t) HN   VN +4 2 ˙ u2 −→ b(t) =  ...  ⇒ ∂u −→ H b(t) ∂t    ∂ 4 u2 VN +4 bN (t) −→ M 4 Hb(t) ∂x4   c0 (t) ∂v1 VN VN   ˙ v1 −→ c(t) =  ...  ⇒ −→ c(t) ∂t cN (t)   d0 (t) ∂v2 VN ˙ VN   −→ d(t) v2 −→ d(t) =  ...  ⇒ ∂t dN (t)   k1,0 VN +4   .. k(x) −→ k1 =   . k1,N +4   k2,0 VN   k(x) −→ k2 =  ...  k2,N 

From (3.3a), we get

T ˙ [Ga(t)] TN +4 = [c(t)]T TN

After taking inner product with Tj (x), j = 0, . . . , N , we get ˙ = cj (t), j = 0, . . . , N [G]j a(t) where cj (t) denotes the j th component of c(t), [G]j is the j th row of G. Writing them in matrix-vector form yields ˙ G1 a(t) = IN +1 c(t) (3.37) where G1 , of order (N + 1) × (N + 1), is a submatrix obtained from the first N + 1 rows of G, and IN +1 is identity matrix of order N + 1. Similarly from (3.3c), it can be shown that ˙ H1 b(t) = IN +1 d(t)

(3.38)

where H1 , of order (N + 1) × (N + 1), is a submatrix obtained from the first N + 1 rows of H. The remaining equations to be discretized are equations (3.3b) and (3.3d). Using (2.4) for functions expanded in the function spaces, VN and VN +4 we get that £ ¤T k(x)u1 = [Wk1 Ga(t)]T TN +4 + Wk∗1 Ga(t) TN∗ +4 £ ¤T k(x)u2 = [Wk1 Hb(t)]T TN +4 + Wk∗1 Hb(t) TN∗ +4 £ ¤T k(x)v1 = [Wk2 c(t)]T TN + Wk∗2 c(t) TN∗ £ ¤T k(x)v2 = [Wk2 d(t)]T TN + Wk∗2 d(t) TN∗ 17

With the function space expansion of u1 , u2 and v1 , we then obtain from (3.3b) that £ ¤T ˙ T TN = −p1 M 4 Ga(t) TN +4 − [Wk1 Ga(t)]T TN +4 [c(t)] £ ¤T − Wk∗1 Ga(t) TN∗ +4 + [Wk1 Hb(t)]T TN +4 £ ¤T + Wk∗1 Hb(t) TN∗ +4 − δ [Wk2 c(t)]T TN £ ¤T −δ Wk∗2 c(t) TN∗ Upon taking inner product with Tj (x), j = 0, . . . , N , we get £ ¤ c˙j (t) = −p1 M 4 G j a(t) − [Wk1 G]j a(t) + [Wk1 H]j b(t) −δ [Wk2 ]j c(t),

j = 0, . . . , N

where [ ∗ ]j denotes the j th row of the indicated matrix, cj (t) denotes the j th component of c(t), dj (t) denotes the j th component of d(t). They can be written in matrix-vector form as, ˙ = A1 a(t) − B1 a(t) + B2 b(t) − δC1 c(t) IN +1 c(t)

(3.39)

where A1 , of order (N + 1) × (N + 1), is a submatrix obtained from the first N + 1 rows of −p1 M 4 G, and B1 , of order (N + 1) × (N + 1), is a submatrix obtained from the first N + 1 rows of Wk1 G, B2 , of order (N + 1) × (N + 1), is a submatrix obtained from the first N + 1 rows of Wk1 H, and C1 , of order (N + 1) × (N + 1), is equal to Wk2 . Similarly with the function space expansion of u1 , u2 and v2 , we then obtain from (3.3d) that, £ ¤ d˙j (t) = −p2 M 4 H j b(t) + [Wk1 G]j a(t) − [Wk1 H]j b(t) −δ [Wk2 ]j d(t),

j = 0, . . . , N

Writing them in matrix-vector form yields, ˙ = A2 b(t) + B1 a(t) − B2 b(t) − δC1 d(t) IN +1 d(t)

(3.40)

where A2 , of order (N + 1) × (N + 1), is a submatrix obtained from the first N + 1 rows of −p2 M 4 H. If we now let   a(t)  b(t)   z(t) =   c(t)  d(t) then (3.37)-(3.40) may be combined to give     G1 0 0 0 0 0 IN +1 0  0 H1  0 0  0 0 0 IN +1   z(t)  z(t)   ˙ =  0  A 1 − B1 0 IN +1 0  B2 −δC1 0  0 0 0 IN +1 B1 A 2 − B2 0 −δC1 ˙ = Dz(t) ≡ Qz(t) (3.41) 18

The above equations are the discretized system of equations associated with the new basis technique while v1 and v2 not satisfying the boundary conditions. The ˆ where eigenvalue problem for these equations is obtained by setting z(t) = eλt z,   ˆ a ˆ   b  ˆ= zˆ =    cˆ  , a dˆ 

  a0  ..  , b .  ˆ= aN

  b0 ..  , cˆ =   .  bN

  c0 ..  , dˆ =   .  cN

 d0 ..  .  dN

therefore we obtain

Qλeλt zˆ = Deλt zˆ ⇒ Q−1 D zˆ = λzˆ

(3.42)

which is the eigenvalue problem corresponding to equation (3.41).

4

Numerical result

All the numerical experiments for equations (3.30), (3.36) and (3.42) are done with Matlab 7.9.0.529(64-bit) on a PC with Windows 7, and Q−1 D is computed through Q \ D, which, according to the Matlab documentations, provides better accuracy than Q−1 D.

4.1

Separated problem

Consider the simplest case of p1 = p2 = 1, k(x) = 0 and δ = 0, then (3.3a)-(3.3d) become ∂u1 = v1 ∂t ∂v1 ∂ 4 u1 = − 4 ∂t ∂x

(4.1) (4.2)

with boundary conditions (3.2a)-(3.2b), and ∂u2 = v2 ∂t ∂v2 ∂ 4 u2 = − 4 ∂t ∂x

(4.3) (4.4)

with boundary conditions (3.2c)-(3.2d). They are actually the same equations but with different boundary conditions, and analytically solvable, it can be shown that they have exactly the same solution except that the boundary conditions (3.2a)-(3.2b) cause to 4 zero eigenvalues. The first non-zero eigenvalue computed with Mathematica 7.0.0(64-bit) is exactly 0 + 5.5933213620153309807i. To do the computation, notice that k(x) = 0 and δ = 0, after plugging them into (3.30), (3.36) and (3.42), B in (3.30) becomes 0, B1 and B2 in (3.36) become 0, B1 19

and B2 in (3.42) become 0, therefore each of (3.30), (3.36) and (3.42) can be divided for the corresponding separated problem. For example, in (3.42), we define ¸ ¸ · · IN +1 0, IN +1 0 , Q2 = Q1 = 0 H1 0 G1 ¸ ¸ · · 0 IN +1 0 IN +1 , D2 = D1 = A2 0 A1 0 The eigenvalue problem associated with (4.1)-(4.2) with boundary conditions (3.2a)(3.2b) is Q1 −1 D1 zˆ = λzˆ and the eigenvalue problem associated with (4.3)-(4.4) with boundary conditions (3.2c)-(3.2d) is Q2 −1 D2 zˆ = λzˆ Table (4.1), (4.2) and (4.3) list the first non-zero eigenvalue given by (3.30), (3.36) and (3.42) when using boundary condition (3.2a)-(3.2b), respectively. Table (4.4), (4.5) and (4.6) list the first non-zero eigenvalue given by (3.30), (3.36) and (3.42) when using boundary condition (3.2c)-(3.2d), respectively. As shown in the tables (4.1)-(4.3), the significant digits in imaginary part of the result given by (3.30) get lost rapidly, but (3.36) and (3.42) can always give at least 10 and 12 significant digits, respectively. As to the real part, the exact value should be 0, hereby (3.42) gives the best result, while (3.36) worse and (3.30) worst. From tables (4.4)-(4.6), we can see that the significant digits in the imaginary part of the result given by (3.30) get lost very rapidly, while (3.36) and (3.42) can always give at least 7 and 9 significant digits, respectively. As to the real part, (3.42) gives the best result, while (3.36) worse and (3.30) worst.

Table 4.1: Separated: result of (3.30) with (3.2a)-(3.2b) # of terms first non-zero eigenvalue 20 -8.865743425859E-014 + 5.593321362019 i 40 3.021959310394E-012 + 5.593321362566 i 60 -4.140124054946E-012 + 5.593321368810 i 80 5.346158007766E-012 + 5.593321352681 i 100 1.437523170803E-011 + 5.593321464868 i 120 4.602572400961E-011 + 5.593321132434 i 140 3.064596370531E-010 + 5.593322449676 i 160 1.671746862752E-010 + 5.593317063761 i 180 -1.924573923972E-010 + 5.593321985781 i 200 9.034823417390E-010 + 5.593315177083 i

20

Table 4.2: Separated: result of (3.36) with (3.2a)-(3.2b) # of terms first non-zero eigenvalue 20 -4.701794509288E-014 + 5.593321362015 i 40 5.602254771198E-013 + 5.593321362016 i 60 1.032139651524E-012 + 5.593321362015 i 80 9.888312391126E-012 + 5.593321362003 i 100 3.653119473590E-012 + 5.593321361995 i 120 1.444921959859E-011 + 5.593321362028 i 140 9.900968933607E-012 + 5.593321362045 i 160 1.439803831715E-011 + 5.593321362063 i 180 1.945150984728E-011 + 5.593321362153 i 200 7.179307148775E-011 + 5.593321362106 i

Table 4.3: Separated: result of (3.42) with (3.2a)-(3.2b) # of terms first non-zero eigenvalue 20 3.157196726278E-016 + 5.593321362016 i 40 -7.251144129583E-014 + 5.593321362015 i 60 4.255485395489E-014 + 5.593321362016 i 80 -4.264470720994E-014 + 5.593321362015 i 100 1.506850200172E-013 + 5.593321362015 i 120 1.480550572549E-014 + 5.593321362015 i 140 -7.669451011771E-014 + 5.593321362015 i 160 -4.677657458423E-014 + 5.593321362015 i 180 -1.236031862663E-013 + 5.593321362016 i 200 -1.205979760499E-014 + 5.593321362014 i

Table 4.4: Separated: result of (3.30) with (3.2c)-(3.2d) # of terms first non-zero eigenvalue 20 2.526550228776E-013 + 5.593321366400 i 40 -4.522384629820E-011 + 5.593319852832 i 60 9.841785656211E-013 + 5.593280899477 i 80 -1.531927046838E-009 + 5.593640098176 i 100 2.496032726277E-010 + 5.596699342572 i 120 9.702530544133E-011 + 5.603683972977 i 140 1.064800329789E-008 + 5.570719168921 i 160 1.623406346150E-008 + 5.780657216792 i 180 -1.490116118673E-008 + 5.675292828342 i 200 -1.239398311875E-008 + 5.665632837290 i

21

Table 4.5: Separated: result of (3.36) with (3.2c)-(3.2d) # of terms first non-zero eigenvalue 20 7.742938235022E-015 + 5.593321362022 i 40 1.283478531788E-014 + 5.593321361845 i 60 4.075147256319E-012 + 5.593321360575 i 80 7.546856477421E-013 + 5.593321361450 i 100 3.710117852046E-012 + 5.593321368787 i 120 1.153386862743E-011 + 5.593321385432 i 140 -1.020835438479E-011 + 5.593321329017 i 160 -1.406817569475E-011 + 5.593321370792 i 180 1.150672454184E-011 + 5.593321430929 i 200 -5.847618981240E-011 + 5.593321475230 i

Table 4.6: Separated: result of (3.42) with (3.2c)-(3.2d) # of terms first non-zero eigenvalue 20 0.000000000000E-000 + 5.593321362014 i 40 -3.928140703880E-016 + 5.593321362005 i 60 8.398799190894E-016 + 5.593321361949 i 80 -1.254545241563E-014 + 5.593321362040 i 100 -1.526556658860E-015 + 5.593321361852 i 120 -3.284910984254E-015 + 5.593321361727 i 140 -1.376676550535E-013 + 5.593321361747 i 160 1.016147140932E-014 + 5.593321361707 i 180 -1.488707161018E-014 + 5.593321362184 i 200 -1.162047888470E-015 + 5.593321361964 i

4.2

Coupled problem 1 2

1 2

Choose parameters k(x) = 0.1(e−16(x− 2 ) + e−16(x+ 2 ) ), p1 = p2 = 1, δ = 1. Table (4.7), (4.8) and (4.9) list the eigenvalue with imaginary part closest to 5.5933213620153309807i given by (3.30), (3.36) and (3.42), respectively. Just as the separated problem, (3.30) loses its significant digits quickly in real part, and gives only 3 significant digits in imaginary part, while (3.36) gives 2 and 6 significant digits in real and imaginary parts, and (3.42) gives 5 and 9 significant digits in real and imaginary parts. Notice that (3.42) is deduced with v1 and v2 not satisfying the boundary conditions, and the numerical results are consistent with this, but it can provide better approximation to the eigenvalues. This is an interesting phenomenon. Figures (4.1)-(4.3) show the distribution of the eigenvalues given by (3.36) with parameters presented in section 4.2. Obviously, all the eigenvalues are in the left half plane, specifically Re(λ) < −0.5E − 4, which indicates the stability of the system simulated by equation (3.36). It is easy to understand, since this is a conservative 22

system with internal energy dissipation, so it must be stable.

Table 4.7: Coupled: result of (3.30) # of terms eigenvalue closest to 5.5933213620153309807i 20 -1.696600794218E-003 + 5.593587490912 i 40 -2.014203095153E-003 + 5.593647464088 i 60 -2.013375812398E-003 + 5.593647220497 i 80 -2.016298462637E-003 + 5.593645643022 i 100 -2.000675748912E-003 + 5.593645196102 i 120 -2.022480683191E-003 + 5.593662912802 i 140 -1.654820503910E-003 + 5.593318134435 i 160 -2.224705539659E-003 + 5.592515222158 i 180 3.560195441660E-003 + 5.582907446855 i 200 -7.681585147621E-004 + 5.591396320852 i

Table 4.8: Coupled: result of (3.36) # of terms eigenvalue closest to 5.5933213620153309807i 20 -3.074759362660E-004 + 5.593372889475 i 40 -3.654157349264E-004 + 5.593383211466 i 60 -3.652607020196E-004 + 5.593383186288 i 80 -3.652532971138E-004 + 5.593383184428 i 100 -3.652398643337E-004 + 5.593383186897 i 120 -3.653078373088E-004 + 5.593383157640 i 140 -3.652671952068E-004 + 5.593383142103 i 160 -3.654157221780E-004 + 5.593383161375 i 180 -3.651049628948E-004 + 5.593383093593 i 200 -3.651957999956E-004 + 5.593384157939 i

23

Table 4.9: Coupled: result of (3.42) # of terms eigenvalue closest to 5.5933213620153309807i 20 -5.580475752882E-004 + 5.593372228902 i 40 -3.653805144792E-004 + 5.593383209869 i 60 -3.652897489307E-004 + 5.593383185712 i 80 -3.652536956824E-004 + 5.593383184645 i 100 -3.652542187657E-004 + 5.593383184697 i 120 -3.652541558668E-004 + 5.593383184712 i 140 -3.652541474563E-004 + 5.593383184708 i 160 -3.652551354990E-004 + 5.593383184569 i 180 -3.652535417538E-004 + 5.593383184839 i 200 -3.652578033432E-004 + 5.593383184065 i

8

3

x 10

2

1

0

−1

−2

−3 −0.012

−0.01

−0.008

−0.006

−0.004

−0.002

0

Figure 4.1: Eigenvalues given by (3.36) with N = 200

4.3

Conclusion

From the comparison of the results for the separated problem, we can see that the new basis technique is more reliable than the tau method, since it provides better numerical stability and better accuracy, thus we think the result given by new basis technique for the coupled problem should also be more promising. This does not mean Tau method is worse, actually when we computed the OrrSommerfeld stability equation posed in reference [2] with new basis technique, we get almost exactly the same result as Tau method. But for the equations presented in this paper, new basis technique give more promising result. 24

4

x 10

3

2

1

0

−1

−2

−3

−5.5396

−5.5394

−5.5392

−5.539

−5.5388

−5.5386

−5.5384

−5.5382 −3

x 10

Figure 4.2: Zoom-in of part of figure (4.1)

6

x 10

8

6

4

2

0

−2

−4

−6

−8

−3.5

−3

−2.5

−2

−1.5

−1

−0.5

0 −4

x 10

Figure 4.3: Zoom-in of part of figure (4.1)

25

This paper demonstrates how to apply the new basis technique to a one-dimension problem, it is valuable to study whether or not it can be applied to higher dimension problem and how to do that.

A

Appendix

Here are some properties of Chebyshev polynomials: T0 (x) T1 (x) Tn (x) Tn (±1)

= = = =

1 x 2xTn−1 (x) − Tn−2 (x), n ≥ 2 (±1)n , n = 0, 1, . . .

dk Tn (±1) = (±1)n+k dxk Particularly,

Qk−1

− m2 ) , k ≥ 1, n = 0, 1, . . . (2k − 1)!!

m=0 (n

2

(A.1)

2 2 2 d3 n+3 n (n − 1)(n − 4) T (±1) = (±1) n dx3 15 2 2 d2 n (n − 1) Tn (±1) = (±1)n+2 2 dx 3 d Tn (±1) = (±1)n+1 n2 dx

Define inner product hu, vi = then

Suppose

Z

1

−1

√

uv dx, 1 − x2

  π, i = j = 0 π hTi (x), Tj (x)i = , i = j 6= 0 2   0, i 6= j f (x) =

N X

N

an Tn (x),

n=0

(k)

X dk f (x) = a(k) n Tn (x) dxk n=0

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 where

( 2, j = 0 cj = 1, j = 1, 2, . . .

26

(A.2)

By (A.2) and mathematical induction, it can be  (k)   a0  a(k)    1  k  . =M   ..   (k)

aN

shown that  a0 a1   ..  .  aN

(A.3)

where matrix M , depending only on N , is generated from (A.2) as below. (1) Consider k = 1. Obviously aN = 0, because the order of the polynomial is decreased by one after taking its first derivative. Therefore the coefficient of the N th order term must be 0. (1) Second, aN −1 = 2N aN , because T0 (x) T1 (x) T2 (x) T3 (x) .. . Tn−1 (x) Tn (x) dTn (x) dx

= = = =

1 x = 21−1 x 2x2 − 1 = 22−1 x2 + . . . 4x3 + . . . = 23−1 x3 + . . .

= 2n−2 xn−1 + . . . = 2n−1 xn + . . . = n ∗ 2n−1 xn−1 + . . . = 2n ∗ 2n−2 xn−1 + . . . = 2nTn−1 (x) + . . .

After taking the first derivative of a polynomial of order N , it becomes a polynomial of order N − 1. From the above calculation the term of order N − 1 is TN −1 (x) which is obtained from the first derivative of aN TN (x). Therefore using the result above, we (1) get that aN −1 = 2N aN . (1) (1) Now that aN = 0 and aN −1 = 2N aN , by (A.2), all the coefficients can be com(1) (1) (1) puted one by one, from aN −2 , aN −3 , . . ., until a0 . P (1) Suppose an = N i=0 sn,i ai , n = 0, . . . , N , then we have 

where

   

(1)

a0 (1) a1 .. .

(1)

aN



  M = 





    =M  

s0,0 s1,0 .. . sN,0

a0 a1 .. .

aN

    

 s0,1 . . . s0,N s1,1 . . . s1,N     sN,1 . . . sN,N 27

B

Matlab Code

B.1

report Listing B.1: report.m

1

% main f i l e

2 3 4

clear ; close a l l ;

5 6 7 8

% N i s b e t we e n n1 and n2 n1 = 2 0 ; n2 = 2 0 0 ;

9 10 11 12

% 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;

13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34

% Tau method : s e p a r a t e d problem f or p = 1 : 2 %boundary c o n d i t i o n : 1 or 2 f or 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 or m = n1 +1: n2 tmp = Tau1Beam (m, r , p ) ; pL = length ( r e s (m−n1 , : ) ) ; f or 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 or m = 1 : pL DrawSave1 ( 0 , r e s ( : ,m) , r , p , 0 , n1 , n2 ) ; end end end

35 36 37 38 39 40 41 42 43 44 45 46 47 48

% Tau method : c o u p l e d problem f or 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 or m = n1 +1: n2 tmp = Tau2Beam (m, r ) ; pL = length ( r e s (m−n1 , : ) ) ; f or 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 ) ) ;

28

end end f or m = 1 : pL DrawSave2 ( 0 , r e s ( : ,m) , r , 1 , n1 , n2 ) ; end

49 50 51 52 53 54

end

55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76

% New b a s i s t e c h : v1 & v2 use boundary c o n d i t i o n : s e p a r a t e d problem f or p = 1 : 2 %boundary c o n d i t i o n : 1 or 2 f or r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = NBUseBC1Beam( 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 or m = n1 +1: n2 tmp = NBUseBC1Beam(m, r , p ) ; pL = length ( r e s (m−n1 , : ) ) ; f or 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 or m = 1 : pL DrawSave1 ( 1 , r e s ( : ,m) , r , p , 0 , n1 , n2 ) ; end end end

77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96

% New b a s i s t e c h : v1 & v2 use boundary c o n d i t i o n : c o u p l e d problem f or r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = NBUseBC2Beam( 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 or m = n1 +1: n2 tmp = NBUseBC2Beam(m, r ) ; pL = length ( r e s (m−n1 , : ) ) ; f or 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 or m = 1 : pL DrawSave2 ( 1 , r e s ( : ,m) , r , 1 , n1 , n2 ) ; end end

97 98 99 100 101 102

% New b a s i s f or p = 1 : 2 f or r = res tmp

t e c h : v1 & v2 no boundary c o n d i t i o n : s e p a r a t e d problem %boundary c o n d i t i o n : 1 or 2 0 : 1 %damping term : 0 or 1 = []; = NBNoBC1Beam( n1 , r , p ) ;

29

tmp = tmp ( imag ( tmp )>=0) ; tmp = tmp ( abs ( tmp ) < r a d i u s ) ; tmp = sort ( tmp ) ; r e s ( 1 , : ) = tmp ; f or m = n1 +1: n2 tmp = NBNoBC1Beam(m, r , p ) ; pL = length ( r e s (m−n1 , : ) ) ; f or 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 or m = 1 : pL DrawSave1 ( 2 , r e s ( : ,m) , r , p , 0 , n1 , n2 ) ; end

103 104 105 106 107 108 109 110 111 112 113 114 115 116

end

117

end

118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138

% New b a s i s t e c h : v1 & v2 no boundary c o n d i t i o n : c o u p l e d problem f or r = 0 : 1 %damping term : 0 or 1 res = [ ] ; tmp = NBNoBC2Beam( 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 or m = n1 +1: n2 tmp = NBNoBC2Beam(m, r ) ; pL = length ( r e s (m−n1 , : ) ) ; f or 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 or m = 1 : pL DrawSave2 ( 2 , r e s ( : ,m) , r , 1 , n1 , n2 ) ; end end

B.2

c Listing B.2: c.m

1 2 3 4 5 6 7 8

% c(k) function r = c ( k ) i f ( k==0) r = 2.0; else r = 1.0; end end

B.3

ChebyshevExpansion Listing B.3: ChebyshevExpansion.m

1

% 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

30

2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18

% 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 or 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

B.4

ChebyshevMultiplication Listing B.4: ChebyshevMultiplication.m

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

% 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) ; f or n = 0 : L f or 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 end end

B.5

GetClosestPoint Listing B.5: GetClosestPoint.m

1 2 3 4 5 6 7 8 9 10

% 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