耦合微分的緊緻數值方法解一階KDV方程

國 立 交 通 大 學

應用數學系

碩 士 論 文

耦合微分的緊緻數值方法

解一階 KDV 方程

Coupled Derivatives Compact Schemes

for One-Dimensional KdV Equation

研 究 生：李雅羚

指導老師：賴明治 教授

耦合微分的緊緻數值方法

解一階 KDV 方程

Coupled Derivatives Compact Schemes

for One-Dimensional KdV Equation

研 究 生: 李雅羚 Student：Ya-Ling Li

指導教授: 賴明治 Advisor：Ming-Chih Lai

國 立 交 通 大 學

應 用 數 學 系

碩 士 論 文

Submitted to Department of Applied Mathematics College of Science

National Chiao Tung University in Partial Fulfillment of the Requirements

for the Degree of Master

in

耦合微分的緊緻數值方法解一階 KDV 方程

學生: 李雅羚 指導教授: 賴明治 教授

國立交通大學應用數學系(研究所)碩士班

摘 要

這篇論文主要之目的是使用耦合微分的緊緻數值方法

來解一階KDV方程。首先，我們先回顧一階和二階耦合微分

的緊緻數值方法。接著，我們會學習一階和三階耦合微分的

緊緻數值方法。再來，我們簡要地介紹Runge-Kutta Methods。

最後，我們會給一些例子並且列出數值結果，然後做出結論。

Coupled Derivatives Compact Schemes

for One-Dimensional KdV Equation

student：Ya-Ling Li Advisor：Dr. Ming-Chih Lai

Department (Institute) of Applied Mathematics

National Chiao Tung University

ABSTRACT

The primary objective of this thesis is to use coupled derivatives

compact schemes (CD) for solving one-dimensional KDV equation. First,

we review the coupled first and second derivatives scheme and then we

study the coupled first and third derivatives scheme. Next, we introduce

roughly the Runge-Kutta methods. Finally, we give some examples and

show numerical results, and the conclusion follows.

誌 謝

這篇論文的完成，首先要感謝我的指導老師 賴明治教

授。在這兩年來，老師除了在學問上的指導令我獲益匪淺之

外，其對於研究事物的態度更是讓人如沐春風，謹此致上我

最誠摯的敬意與謝意。口試期間，也承蒙林吉田老師、張書

銘老師及鄭博文老師費心審閱並提供了寶貴的意見，使得本

論文得以更加的完備，永誌於心。

在這兩年求學的過程中，感謝曾昱豪學長在我遇到問題

時，總是細心及耐心地給予我意見及靈感的啟發；感謝吳恭

儉學長和曾孝捷同學在我遇到瓶頸時給予無私的協助與關

心；另外，也要感謝陳冠羽學弟以及實驗室之前的學長姐們

給我的支持與鼓勵。

除此之外，更要感謝我研究室的同學們:倖綺、園芳、

千砡、碧維、慧萍以及怡菁，謝謝她們讓我的研究所生活更

多采多姿，有了她們的陪伴與支持，我才有這些美好的回憶。

最後，要感謝我的家人們，以及男友謝俊鴻的默默關懷

與加油打氣。願與所有關心我的人一起分享這份喜悅，再次

地感謝所有幫助過我及關心過我的人，謝謝你們！

目

錄

中文提要

………

i

英文提要

………

ii

誌謝

………

iii

目錄

………

iv

一、

Introduction ………

1

二、

2.1

2.2

The Coupled First and Second Derivatives Scheme …

First Equation for First and Second Derivatives ……

Second Equation for First and Second Derivatives …

2

3

4

三、

3.1

3.2

3.3

四、

4.1

4.2

4.3

The Coupled First and Third Derivatives Scheme …

First Equation for First and Third Derivatives ……

Third Equation for First and Third Derivatives ……

The Scheme ………

Runge-Kutta Methods ………

Euler’s Method ………

Second-Order Runge-Kutta Method ………

Fourth-Order Runge-Kutta Method ………

6

7

8

8

11

11

11

12

五、

5.1

5.2

5.3

Numerical Examples ………

Example 1 ………

Example 2 ………

Example 3 ………

13

13

14

18

六、

Conclusion ………

21

References

………

22

Coupled Derivatives Compact Schemes for

One-Dimensional KdV Equation

Student : Ya-Ling Li

Advisor : Ming-Chih Lai

1

Introduction

The KdV (Korteweg-de Vries) equation is a nonlinear partial differential equation. Two Dutch mathematicians D. J. Korteweg and G. deVries discovered this famous KdV equa-tion

ut+ uux+ uxxx = 0

when they derived the shallow water wave, and they only considered dispersion but ignored the dissipation of the energy.

The KdV equation appears in a great number of physical situations. The reason for the ubiquitous incident of KdV equation is at least twofold. First, it associates simple dispersion with weak nonlinearity. Second, using the asymptotic method of manifold scales it can be shown that it describes (on suitable scales) the Riemann invariants of any hyperbolic system with weak nonlinearity and dispersion. The KdV equation is integrable [1], i.e. it can be written as the compatibility condition of a pair of linear eigenvalue equations, called the Lax pair [2]. In many physical situations, KdV appears in the form of an initial-boundary value problem on the semi-infinite line. This is for example the case of a certain laboratory study of water waves [3]. However, the solution of initial-boundary value problems for integrable equations was until recently open.

The KdV equation is the original of integrable equations. Initially presented as an equation with a solitary wave-type solution, it turned out (much) later to hold solutions with an arbitrary number of flexibly scattering solitary waves. The latter were observed numerically by Kruskal and Zabusky, who decided to call this type of solitary waves ’soli-tons’. In the years that followed this discovery, Kruskal and his collaborators went on to show that the KdV equation held an infinite number of conservation laws and, as the eventual explanation of these properties, they produced a linear differential system the compatibility of which is just the KdV equation. This linear system is traditionally called the Lax pair and allows the efficient linearization of the nonlinear equation. Using tech-niques developed in the theory of the inverse scattering problems in quantum mechanics (reconstruction of the potential from the scattering data) one can reduce the solution of the KdV equation to that of a linear integrodifferential one. This was the final proof of the integrability of KdV. The discovery of an integrable PDE and the techniques for its integration opened a whole new domain that is still the center of extreme activity a quarter-century later.

This paper presents a family of finite difference schemes for the first and third deriva-tives of smooth functions. The schemes are Hermitian and symmetric. The objective of this paper is to develop this family of schemes and to assess their potential for computa-tions of the KdV equation. The schemes will be referred to as the “coupled-derivative,” or “C-D” schemes.

2

The Coupled First and Second Derivatives Scheme

First, we review a family of finite schemes for the first and second derivatives of smooth functions [5]. The schemes are Hermitian and symmetric. When defined on a uniform mesh, the schemes are of the form

a1f 0 i−1+ a0f 0 i + a2f 0 i+1+ h(b1f 00 i−1+ b0f 00 i + b2f 00 i+1) = 1

h(c1fi−2+ c2fi−1+ c0fi+ c3fi+1+ c4fi+2). (1) Throughout this paper, h denotes the uniform mesh spacing. The interior scheme is of the form given by Eq.(1). Simultaneous solving for f_i0 and f_i00 implies that the number of unknowns is equal to 2M . A total of 2M equations are therefore needed to close the system. Equation(1) may be used to derive two linearly independent equations at each node. This is done as follow. Both sides of Eq.(1) are first expanded in a Taylor series. The resulting coefficients are then matched, such that Eq.(1) maintains a certain order of accuracy. Note that Eq.(1) has 11 coefficients, of which one is arbitrary; i.e., Eq.(1) may be divided through by one of the constants without loss of generality. A convenient choice of the normalization constant is either of a0 or b0. It will be seen that the equation obtained by setting a0 equal to 1 is linearly independent of the equation obtained when b0 is set equal to 1. The two equations may therefore be applied at each node, and the resulting system of 2M equations solved for the nodal values of the first and second derivative. The process of obtaining the two equations is outlined in Sections 2.1 and 2.2.

TABLE I

Taylor Table for a0 = 1 LHS RHS fi 0 c0 f_i0 1+2a1 2(2c4+c3) f_i00 b0 0 f_i000 2h2_(a 1/2! + b2) 2h2(23c4+ c3)/3! f_iiv 0 0 fv i 2h4(a1/4! + b2/3!) 2h4(25c4+ c3)/5! fvi i 0 0 f_ivii 2h6(a1/6! + b2/5!) 2h6(27c4+ c3)/7! fviii i 0 0 fix i 2h8(a1/8! + b2/7!) 2h8(29c4+ c3)/9!

2.1

First Equation(a

= 1) for First and Second Derivatives

Consider first the case where a0 = 1. The symmetry of the schemes requires that a1 = a2, b1 = −b2, c1 = −c4, and c2 = −c3. Equation(1) therefore reduces to

a1f 0 i−1+ f 0 i + a1f 0 i+1+ h(−b2f 00 i−1+ b0f 00 i + b2f 00 i+1) = 1

h[c0fi+ c3(fi+1− fi−1) + c4(fi+2− fi−2)]. (2) Expanding both sides of Eq.(2) in a Taylor series and collecting terms of the same order yields Table I. Note that “LHS” and “RHS” denote the coefficients of fk

i on the left- and right-hand sides, respectively, of Eq.(2).

The Taylor table shows that b0 = c0 = 0. This leaves four undetermined constants(a1, b2, c3, and c4). Expressions for these constants may be obtained by matching the terms in the Taylor table.

When a0 = 1, b0 = 0:

Matching terms up to f_i00 yields

a1 = − 1

2 + c3+ 2c4, b2 arbitrary . Matching terms up to f_iiv yields

a1 = − 1 2 + c3+ 2c4, b2 = 1 12[3 − 4(c3− c4)]. Matching terms up to fvi i yields a1 = 7 16− 15 4 c4, b2 = 1 16(−1 + 36c4), c3 = 15 16− 23 4 c4. Matching terms up to fviii

i yields a1 = 17 36, b2 = − 1 12, c3 = 107 108, c4 = − 1 108. ⇒ 17 36f 0 i−1+ f 0 i + 17 36f 0 i+1+ h( 1 12f 00 i−1− 1 12f 00 i+1) = 1 h[ 107 108(fi+1− fi−1) − 1 108(fi+2− fi−2)]. i.e.

51f_i−10 + 108f_i0+ 51f_i+10 + 9h(f_i−100 − f_i+100 ) = 107

h (fi+1− fi−1) −

fi+2− fi−2 h .

2.2

Second Equation(b

= 1) for First and Second Derivatives

Consider the case where b0 = 1. Note that a tilde is used above the constants to indicate their difference from the constants obtained when a0 = 1; e.g., b1 is replaced by ˜b1. Symmetry requires that ˜b1 = ˜b2, ˜c1 = ˜c4, ˜c2 = ˜c3, and ˜a1 = − ˜a2. Equation(1) therefore becomes ˜ a0f 0 i + ˜a2(f 0 i+1− f 0 i−1) + h( ˜b1f 00 i−1+ f 00 i + ˜b1f 00 i+1) = 1

h[ ˜c1(fi−2+ fi+2) + ˜c2(fi−1+ fi+1) + ˜c0fi]. (3) Expanding both sides of the above equation in a Taylor series and collecting terms of the same order yields the Taylor Table II.

Table II shows that ˜a0 is required to be zero if ˜b0 is equal to one. The resulting equa-tion may therefore be considered an expression for the second derivative. We have five unknown constants( ˜c0, ˜c1, ˜c2, ˜a2, and ˜b1). These constants may be obtained by matching the terms in the above Taylor table and solving the resulting equations.

TABLE II

Taylor Table for b0 = 1 LHS RHS fi 0 c˜0+ 2 ˜c1+ 2 ˜c2 f_i0 a˜0 0 f_i00 h(2 ˜a2 + 2 ˜b1 + 1) 2h(22c˜1+ ˜c2)/2! f_i000 0 0 fiv i 2h3( ˜a2/3! + ˜b1/2!) 2h3(24c˜1+ ˜c2)/4! f_iv 0 0 fvi i 2h5( ˜a2/5! + ˜b1/4!) 2h5(26c˜1+ ˜c2)/6! f_ivii 0 0 f_iviii 2h7( ˜a2/7! + ˜b1/6!) 2h7(28c˜1+ ˜c2)/8! fix i 0 0 f_ix 2h9( ˜a2/9! + ˜b1/8!) 2h9(210c˜1+ ˜c2)/10!

When b0 = 1, a0 = 0:

Matching terms up to f_i00 yields

c0 = −2(c1+ c2), a2 = 1

2(−1 − 2b1+ 4c1+ c2). Matching terms up to f_iiv yields

c0 = −2(c1+ c2), a2 = − 3 4 + c1+ 5 8c2, b1 = 1 4 + c1− c2 8. Matching terms up to fvi i yields c0 = −6 + 54c1, c2 = 3 − 28c1, a2 = 9 8− 33 2 c1, b1 = − 1 8+ 9 2c1. Matching terms up to fviii

i yields c0 = − 13 2 , c1 = − 1 108, c2 = 88 27, a2 = 23 18, b1 = − 1 6. ⇒ 23 18(f 0 i+1−f 0 i−1)+h(− 1 6f 00 i−1+f 00 i − 1 6f 00 i+1) = 1 h[− 1 108(fi−2+fi+2)+ 88 27(fi−1+fi+1)− 13 2 fi]. i.e.

138(f_i+10 − f_i−10 ) − h(18f_i−100 − 108f_i00+ 18f_i+100 ) = −fi+2+ fi−2 h +

352

h (fi+1+ fi−1) − 702

3

The Coupled First and Third Derivatives Scheme

Now, we present a family of finite difference schemes for the first and third derivatives of smooth functions. The schemes are also Hermitian and symmetric. When defined on a uniform mesh, the schemes are of the form

a1f 0 i−1+ a0f 0 i + a2f 0 i+1+ h 2_(b 1f 000 i−1+ b0f 000 i + b2f 000 i+1) = 1

h(c1fi−2+ c2fi−1+ c0fi+ c3fi+1+ c4fi+2). (4) The interior scheme is of the form given by Eq.(4). Simultaneous solving for f_i0 and f_i000 implies that the number of unknowns is equal to 2M . A total of 2M equations are therefore needed to close the system. Equation(4) may be used to derive two linearly independent equations at each node. This is done as follow. Both sides of Eq.(4) are first expanded in a Taylor series. The resulting coefficients are then matched, such that Eq.(4) maintains a certain order of accuracy. Note that Eq.(4) has 11 coefficients, of which one is arbitrary; i.e., Eq(4) may be divided through by one of the constants without loss of generality. A convenient choice of the normalization constant is either of a0 or b0. It will be seen that the equation obtained by setting a0 equal to 1 is linearly independent of the equation obtained when b0 is set equal to 1. The two equations may therefore be applied at each node, and the resulting system of 2M equations solved for the nodal values of the first and third derivative. The process of obtaining the two equations is outlined in Sections 2.3 and 2.4.

TABLE III Taylor Table for a0 = 1

LHS RHS fi 0 c0 f_i0 a0+2a1 2(2c4+c3) f_i00 0 0 f_i000 2h2(a1/2! + b2) 2h2(23c4+ c3)/3! fiv i 0 0 fv i 2h4(a1/4! + b2/2!) 2h4(25c4+ c3)/5! f_ivi 0 0 fvii i 2h6(a1/6! + b2/4!) 2h6(27c4+ c3)/7!

3.1

First Equation(a

= 1) for First and Third Derivatives

Consider first the case where a0 = 1. The symmetry of the schemes requires that a1 = a2, b1 = b2, c1 = −c4, and c2 = −c3. Equation(1) therefore reduces to

a1f 0 i−1+ f 0 i + a1f 0 i+1+ h 2_(b 2f 000 i−1+ b0f 000 i + b2f 000 i+1) = 1

h[c0fi+ c3(fi+1− fi−1) + c4(fi+2− fi−2)]. (5) Expanding both sides of Eq.(5) in a Taylor series and collecting terms of the same order yields Table III. Note that “LHS” and “RHS” denote the coefficients of fk

i on the left-and right-hleft-and sides, respectively, of Eq.(5).

The Taylor table shows that c0 = 0. This leaves four undetermined constants(a1, b2, c3, and c4). Expressions for these constants may be obtained by matching the terms in the Taylor table.

When a0 = 1, b0 = 0:

Matching terms up to f_i0 yields

a1 = 1

2(−1 + 2c3+ 4c4), b2 arbitrary . Matching terms up to f_i000 yields

a1 = 1 2(−1 + 2c3+ 4c4), b2 = 1 12(3 − 4c3+ 4c4). Matching terms up to fv i yields a1 = 1 32(9 + 60c4), b2 = 1 96(−1 + 36c4), c3 = 1 32(25 − 4c4). Matching terms up to f_ivii yields

a1 = 11 48, b2 = −1 48, c3 = 113 144, c4 = −1 36. ⇒ 11 48f 0 i−1+ f 0 i+ 11 48f 0 i+1+ h2( −1 48f 000 i−1− 1 48f 000 i+1) = 1 h( 1 36fi−2− 113 144fi−1+ 113 144fi+1− 1 36fi+2). i.e.

11f_i−10 + 48f_i0+ 11f_i+10 + h2(−f_i−1000 − f_i+1000 ) = 1

3.2

Third Equation(b

= 1) for First and Third Derivatives

When b0 = 1, a0 = 0:

Matching terms up to f_i000 yields

a1 = c3+ 2c4, b2 = 1 6(−3 − 2c3+ 2c4). Matching terms up to fv i yields a1 = −15 8 (1 − c4), b2 = 1 8(1 + 3c4), c3 = −1 8 (15 + c4). Matching terms up to f_ivii yields

a1 = −35 32 , b2 = 9 32, c3 = −185 96 , c4 = 5 12. ⇒ −35 32 f 0 i−1− 35 32f 0 i+1+h 2₍ 9 32f 000 i−1+f 000 i + 9 32f 000 i+1) = 1 h( −5 12fi−2+ 185 96 fi−1− 185 96 fi+1+ 5 12fi+2). i.e.

−35f_i−10 − 35f_i+10 + h2(9f_i−1000 + 32f_i000+ 9f_i+1000 ) = 1

3h(−40fi−2+ 185fi−1− 185fi+1+ 40fi+2).

3.3

The Scheme

The interior scheme involves applying the equations derived in section 2.3 and 2.4 at each node. The resulting system of 2M equations is then solved to obtain f_i0 and f_i000.

Schemes:

11f_i−10 + 48f_i0+ 11f_i+10 + h2(−f_i−1000 − f_i+1000 ) = 1

3h(4fi−2− 113fi−1+ 113fi+1− 4fi+2).

−35f_i−10 − 35f_i+10 + h2(9f_i−1000 + 32f_i000+ 9f_i+1000 ) = 1

3h(−40fi−2+ 185fi−1− 185fi+1+ 40fi+2). Note that the first and third derivatives are coupled in the C-D schemes. The vector of unknowns is therefore equal to [· · · , f_i0, · · · , f_i000, · · · ]T_.

Sixth-Order Scheme:Periodic

The sixth-order scheme on a periodic domain is given by

                              48 11 0 · · · 0 11 | 0 −h2 _{0 · · · · · · 0 −h}2 11 48 11 0 · · · 0 | −h2 _{0 −h}2 _{0 · · · · · · 0} 0 11 48 11 0 · · · 0 | 0 −h2 _{0 −h}2 _{0 · · · 0} .. . . .. . .. . .. ... | ... . .. . .. . .. ... .. . . .. . .. . .. 0 | ... . .. . .. . .. 0 0 . .. . .. . .. 11 | 0 . .. . .. . .. −h2 11 0 · · · 0 11 48 | −h2 _{0 · · · · · · 0 −h}2 ₀ −− −− −− −− −− −− −− −− −− −− −− −− −− −− −− 0 −35 0 · · · 0 −35 | 32h2 _9h2 _{0 · · · · · · 0 9h}2 −35 0 −35 0 · · · 0 | 9h2 32h2 9h2 0 · · · 0 0 0 −35 0 −35 0 · · · 0 | 0 9h2 _32h2 _9h2 _{0 · · · 0} .. . . .. . .. . .. ... | ... . .. . .. . .. ... .. . . .. . .. . .. 0 | ... . .. . .. . .. 0 0 . .. . .. . .. −35 | 0 . .. . .. . .. 9h2 −35 0 · · · 0 −35 0 | 9h2 0 · · · 0 9h2 32h2                                                               f₁0 .. . .. . f_i0 .. . .. . f_M0 −− f₁000 .. . .. . f_i000 .. . .. . f_M000                                 = 1 3h                                 4fM −1− 113fM + 113f2− 4f3 .. . .. .

4fi−2− 113fi−1+ 113fi+1− 4fi+2 .. . .. . 4fM −2− 113fM −1+ 113f1− 4f2 − − − − − − − − − − − − −− −40fM −1+ 185fM − 185f2+ 40f3 .. . .. .

−40fi−2+ 185fi−1− 185fi+1+ 40fi+2 .. . .. . −40fM −2+ 185fM −1− 185f1+ 40f2                                

Sixth-Order Scheme:Nonperiodic

The sixth-order scheme on a nonperiodic domain is given by

                         48 11 | 0 −h2 11 48 11 | −h2 _{0 −h}2 11 48 11 | −h2 _{0 −h}2 . .. . .. . .. | . .. . .. . .. . .. . .. 11 | . .. . .. −h2 11 48 | −h2 ₀ −− −− −− −− −− −− −− −− −− −− −− −− −− 0 −35 | 32h2 _9h2 −35 0 −35 | 9h2 32h2 9h2 −35 0 −35 | 9h2 _32h2 _9h2 . .. . .. . .. | . .. . .. . .. . .. . .. −35 | . .. . .. 9h2 −35 0 | 9h2 _32h2                                                    f₁0 .. . .. . f_i0 .. . f_M0 −− f₁000 .. . .. . f_i000 .. . f_M000                           = 1 3h                           4f−1− 113f0+ 113f2− 4f3 .. . .. .

4fi−2− 113fi−1+ 113fi+1− 4fi+2 .. . 4fM −2− 113fM −1+ 113fM +1− 4fM +2 − − − − − − − − − − − − −− −40f−1+ 185f0− 185f2+ 40f3 .. . .. .

−40fi−2+ 185fi−1− 185fi+1+ 40fi+2 .. . −40fM −2+ 185fM −1− 185fM +1+ 40fM +2                          

4

Runge-Kutta Methods

Let us review the Runge-Kutta methods [13].

4.1

Euler’s Method

The Taylor-series method with n = 1 is called Euler’s method. It looks like this:

x(t + h) = x(t) + hf (t, x).

This formula has the obvious advantage of not requiring any differentiation of f . This advantage is offset by the necessity of taking small values for h to gain acceptable precision. Still, the method serves as a useful example and is of great importance theoretically because existence theorems can be based on it.

4.2

Second-Order Runge-Kutta Method

Let us begin with the Taylor series for x(t + h) :

x(t + h) = x(t) + hx0(t) +h 2 2!x 00 (t) +h 3 3!x 000 (t) + · · · (6)

From the differential equation, we have

x0(t) = f

x00(t) = ft+ fxx0 = ft+ fxf

x000(t) = ftt+ ftxf + (ft+ fxf )fx+ f (fxt+ fxxf ) ..

.

Here subscripts denote partial derivatives, and the chain rule of differentiation is used repeatedly. The first three terms in Equation(6) can be written now in the form

x(t + h) = x + hf + 1 2h 2_(f t+ f fx) + O(h3) = x + 1 2hf + 1 2h[f + hft+ hf fx] + O(h 3_{). (7)}

where x means x(t), f means f (t, x), and so on. We are able to eliminate the partial derivatives with the aid of the first few terms in the Taylor series in two variables:

f (t + h, x + hf ) = f + hft+ hf fx+ O(h2).

Equation(7) can be rewritten as

x(t + h) = x + 1 2hf +

1

2hf (t + h, x + hf ) + O(h 3_).

Hence, the formula for advancing the solution is

x(t + h) = x(t) + h

2f (t, x) + h

or equivalently, x(t + h) = x(t) +1 2(F1+ F2) (8) where  F1 = hf (t, x) F2 = hf (t + h, x + F1)

This formula can be used repeatedly to advance the solution one step at a time. It is called a second-order Runge-Kutta method. It is also known as Heun’s method.

In general, second-order Runge-Kutta formulas are of the form

x(t + h) = x + w1hf + w2hf (t + αh, x + βhf ) + O(h3), (9) where w1, w2, α, and β are parameters at our disposal. Equation(9) can be rewritten with the aid of the Taylor series in two variables as

x(t + h) = x + w1hf + w2h[f + αhft+ βhf fx] + O(h3). (10) Comparing Equations(7) and (10), we see that we should impose these conditions:

   w1+ w2 = 1 w2α = 1₂ w2β = 1₂ (11)

One solution is w1 = w2 = 1₂, α = β = 1, which is the one corresponding to Heun’s method in Equation(8). The system of Equation(11) has solutions other than this one, such as the one obtained by letting w1 = 0, w2 = 1, α = β = 1₂. The resulting formula from(9) is called the modified Euler method :

x(t + h) = x(t) + F2, where

 F1 = hf (t, x)

F2 = hf (t + 1₂h, x + 1₂F1)

Compare this to the standard Euler method, described in Section 3.1.

4.3

Fourth-Order Runge-Kutta Method

The higher-order Runge-Kutta formulas are very tedious to derive, and we shall not do so. The formulas are rather elegant, however, and are easily programmed once they have been derived. Here are the formulas for the classical fourth-order Runge-Kutta method :

5

The Numerical Examples

5.1

Example 1

Consider the equation,

u = cos 4x, x ∈ [0, 2π].

To compute u0 & u000 and compare the absolutely errors with the exact ones,  u0 = −4 sin 4x

u000 = 64 sin 4x

Mesh cond(A) Maxnorm-Error order 16 215.965244021860 3.442540136851857E-003 -32 799.575261801726 7.867157863383767E-006 8.7734 64 3134.01533292119 2.683871747066746E-008 8.1954 128 12471.7756173990 1.013642503266965E-010 8.0486

Table IV : The error for u0.

Mesh cond(A) Maxnorm-Error order 16 215.965244021860 0.356972051819817 -32 799.575261801726 4.292639448863156E-003 6.3778 64 3134.01533292119 6.327597368027682E-005 6.0841 128 12471.7756173990 9.748273868126489E-007 6.0204

Table V : The error for u000.

Mesh cond(A) cond(A)/(Mesh2) 16 215.965244021860 0.843614234460390625 32 799.575261801726 0.780835216603248046875 64 3134.01533292119 0.76514046213896240234375 128 12471.7756173990 0.76121677352288818359375

Table VI : The relation between cond(A) and M esh.

From the Table III, we can know that the order of the scheme is equal to ninth order. So the order of the first derivative should be eighth order, and the order of the third derivative should be sixth order.

Note that the condition number is large, so A is ill conditioned and any numerical solution of Ax = b must be accepted with a great deal of skepticism.

5.2

Example 2

Consider the linear equation,

 ut+ ux+ uxxx = 0

u(x, 0) = sin 2x, x ∈ [0, 2π]

with the exact solution of uexact(x, t) = sin(2(x + 3t)) when T = 1.

0 1 2 3 4 5 6 7 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Mesh u & u exact u uexact

Figure 1: u and uexact for RK3 (when M=16).

0 1 2 3 4 5 6 7 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Mesh u & u exact u u_exact

0 1 2 3 4 5 6 7 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Mesh u & u exact u u_exact

Figure 3: u and uexact for RK4 (when M=16).

0 1 2 3 4 5 6 7 −1 −0.8 −0.6 −0.4 −0.2 0 0.2 0.4 0.6 0.8 1 Mesh u & u exact u u_exact

Mesh Time step Maxnorm-Error order 16 132 unstable -32 1056 unstable -64 8454 unstable -128 67636 unstable

-Table VII : The error for RK1(When T=1 (∆t ∼ ∆x3_)).

Mesh Time step Maxnorm-Error order 16 132 unstable -32 1056 unstable -64 8454 unstable -128 67636 unstable

-Table VIII : The error for RK2(When T=1 (∆t ∼ ∆x3_)).

Mesh Time step Maxnorm-Error order 16 132 5.028345551417179E-004 -32 1056 7.851200683139936E-006 6.0010 64 8454 1.213175020303714E-007 6.0161 128 67636 1.896873163750867E-009 5.9990

Table IX : The error for RK3(When T=1 (∆t ∼ ∆x3)).

Mesh Time step Maxnorm-Error order 16 132 5.104998880844369E-004 -32 1056 7.846674022080058E-006 6.0237 64 8454 1.213253499471323E-007 6.0151 128 67636 1.896881306542864E-009 5.9991

Table X : The error for RK4(When T=1 (∆t ∼ ∆x3_)).

In order to be stable, it is necessary to assume that µ = _(∆x)∆t3 ≤

1

2. So we take ∆t ≈ (∆x)3

8 . Note that T = n · ∆t, where n is the time step, so n depends on ∆t when we set T = 1. From the above tables, we know that it is about sixth order. This is because it is time dependence, so the order of the error depends on the order of the third derivative.

Mesh Time step Maxnorm-Error order 16 10000 0.1835957083217821E-02 -32 10000 unstable -64 10000 unstable -128 10000 unstable

-Table XI : The error for RK1(When T=1 (∆t = 0.0001)).

Mesh Time step Maxnorm-Error order 16 10000 0.5110460338352829E-03 -32 10000 0.7490918246130795E-05 6.0922 64 10000 unstable -128 10000 unstable

-Table XII : The error for RK2(When T=1 (∆t = 0.0001)).

Mesh Time step Maxnorm-Error order 16 10000 0.5113914889436910E-03 -32 10000 0.7848784020964006E-05 6.0258 64 10000 0.1213325657306585E-06 6.0154 128 10000 unstable

-Table XIII : The error for RK3(When T=1 (∆t = 0.0001)).

Mesh Time step Maxnorm-Error order 16 10000 0.5113915040849681E-03 -32 10000 0.7848778151214875E-05 6.0258 64 10000 0.1213372820135783E-06 6.0154 128 10000 unstable

-Table XIV : The error for RK4(When T=1 (∆t = 0.0001)).

From the previous page, we know that ∆t ∼ (∆x)3_{. Besides, since ∆x =} 2π

M, so ∆x ∼ 1 M, then we can get ∆t ∼ (_M1)3_{. Since ∆t = 0.0001, so when M is getting larger, the error} becomes unstable. From these tables above, we can conclude that RK4 is the most stable method and its solution is the most accurate one.

5.3

Example 3

Consider the non-linear equation,

 ut+ u · ux+ uxxx = 0

u(x, 0) = 3 · (sech(x/2))2, x ∈ [0, 2π]

with the exact solution of uexact(x, t) = 3 · (sech((x − t)/2))2 when T = 1.

0 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 Mesh u & u exact u u_exact

Figure 5: u and uexact for RK3 (when M=16).

0 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 Mesh u & u exact u u_exact

0 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 Mesh u & u exact u u_exact

Figure 7: u and uexact for RK4 (when M=16).

0 1 2 3 4 5 6 7 0 0.5 1 1.5 2 2.5 3 Mesh u & u exact u u_exact

Mesh Time step Maxnorm-Error order 16 158 unstable -32 1159 unstable -64 8857 unstable -128 69234 unstable

-Table XV : The error for RK1(When T=1 (∆t ∼ ∆x3_)).

Mesh Time step Maxnorm-Error order 16 158 unstable -32 1159 unstable -64 8857 unstable -128 69234 unstable

-Table XVI : The error for RK2(When T=1 (∆t ∼ ∆x3_)).

Mesh Time step Maxnorm-Error order 16 158 2.583180630577608E-005 -32 1159 5.098193027186504E-007 5.6630 64 8857 8.033804399509847E-009 5.9878 128 69234 1.079924483171624E-010 6.2171

Table XVII : The error for RK3(When T=1 (∆t ∼ ∆x3)).

Mesh Time step Maxnorm-Error order 16 158 2.524653262514498E-005 -32 1159 5.132999952306427E-007 5.6201 64 8857 8.047767008356743E-009 5.9951 128 69234 1.079270006698607E-010 6.2205

Table XVIII : The error for RK4(When T=1 (∆t ∼ ∆x3_)).

Observe this non-linear case, we can get that it is about sixth order. This is because it is time dependence, so the order of the error depends on the order of the third derivative.

6

Conclusions

A family of finite difference schemes for the first and third derivatives of smooth functions were derived. We have extended it to the KdV equation. The schemes are Hermitian and symmetric. They are different from the schemes in that the first and third derivatives are simultaneously evaluated.

Consider that the KdV equation requires both first and third derivatives of the vari-ables, the proposed schemes appear to be attractive alternatives to the schemes which the first and third derivatives are simultaneously evaluated for computations of the KdV equation.

References

[1] C. S. Gardner, J. M. Green, M. D. Kruskal, and R. M. Miura: Phys. Rev. Lett. v.19, pp. 1095(1967). For recent developments, see A. S. Fokas and V. E. Za-kharov(eds),Important Developments in Soliton Theory, Springer-Verlag, New York, 1993.

[2] P, Lax.: Comm. Pure Appl. Math. v.21, pp. 467(1968).

[3] J. Bona, W. G. Pritchard, and L. R. Scott: Philos. Trans. Royal Soc. London, A v.302, pp. 457-510(1981).

[4] S.K. Lele: Compact finite difference schemes with spectral-like resolution, J.Comput.Phys. v.103, pp. 16(1992).

[5] K. Mahesh: A family of high order finite difference schemes with good spectral resolution, J.Comput.Phys. v.145, pp. 332(1998).

[6] M. Ben-Artzi., J.-P. Croisille, D. Fishelov, and S. Trachtenberg: A pure-compact scheme for the streamfunction formulation of Navier-Stokes equations, J.Comput.Phys. v.205, pp. 640(2005).

[7] T. Nihei and K. Ishii: A fast solver of the shallow water equations on a sphere using a combined compact difference scheme, J.Comput.Phys. v.187, pp. 639(2003).

[8] P.C. Chu and C. Fan: A three-point combined compact difference scheme, J.Comput.Phys. v.140, pp. 370(1998).

[9] W. E and J.-G. Liu: Essentially compact schemes for unsteady viscous incompressible flows, J.Comput.Phys. v.126, pp. 122(1996).

[10] A.S. Fokas: The Korteweg-da Vries Equation and Beyond, Acta Applicandae Math-ematicae v.39 pp. 295-305(1995).

[11] B. Grammaticos and V. Papageorgiou: KdV Equations and Integrability Detectors, Acta Applicandae Mathematicae v.39, pp. 335-348(1995).

[12] Lloyd N. Trefethen: Spectral Methods in MATLAB, pp. 108-111.

[13] D. Kincaid and W. Cheney: Numerical Analysis:Mathematics of Scientific Comput-ing, Third Edition(2002).