One Dimensional Heat Equation

ABSTRACT

CHAPTER ONE

INTRODUCTION

  It’s deeply truth that the search for the exact solution in our world-problems is needed for each of us, but unfortunately, not all problems can be solved exactly; because of nonlinearity and complicated geometry. In many interesting applications, exact solutions may be impossible and unattainable, so we need to approximate. One of the most beautiful branch in mathematics science can deal with these kind of problems, which is (numerical methods). From the earliest times, the development of the high speed-digital computers and modern electronic devices is widely enhances the use of numerical methods in many various of branch in science and engineering (Fausett, 2002; Ozisik, 1994).

  In this dissertation we study the application (in particular accuracy and convergence) of two kinds of numerical methods based on Adomian Decomposition Method (ADM) and Forward Time Central Space (FTCS) method when applied to the heat equation.

Heat equation is an important partial differential equation (pde) used to describe various phenomena in many applications of our daily life.

1.1 Numerical methods

  One of the earliest mathematical writers in this field was by the Babylonians (3,700 years ago). There is evidence that they knew how to find the numerical solutions for quadratic equations, also they approximated the root for an integer number. Numerical methods are: the methods often used to solve ordinary differential equations (ode’s) and partial differential equations (ode’s) (Fauset, 2002).

  The finite difference method and finite element method are an example of the numerical methods.

The finite difference method will be considered in more details in chapter 2.

PDE’s (background)

 Many physical phenomena can be described by partial differential equations (pde’s), to better understanding the phenomena, we need to formulate and solve it using the pde’s. so , what is the pde’s?

Partial differential equations (PDE’s) : are the class of differential equations involving unknown functions of several variables (Wong Y.Y,W,.T.C,J.M,2005).

A pde for a function u(x,y) can be written as:

The order of the pde is the highest order of its derivatives.

  There are three types of the pde’s, elliptic, parabolic and hyperbolic, where they are classified in order of their discriminates.

For example, consider the second order pde in 2 variables:

Auxx+Buxy+Cuyy=f(x,y,u,ux,uy)

Where A, B, C and f are functions of u. So the pde is said to be:

Hyperbolic pde if: B2-4AC>0 .For example:A  uxx-utt=0 .(wave equation) .

Elliptic pde if : B2-4AC<0 .For example uxx+utt=0.(laplace equation)

Parabolic pde if : B2-4AC=0.For example uxx-ut=0.(Heat equation).  (Wong Y.Y,W,.T.C,J.M,2005).

Heat Equation

  One dimensional heat equation that will be used in this dissertation is:

Where is positive constant. In order for this equation to be solved, the initial conditions (IC) and the boundary conditions (BC) should be found. (Wong Y.Y,W,.T.C,J.M,2005).

  The one dimensional heat equation describes the distribution of heat, heat equation almost known as diffusion equation; it can arise in many fields and situations such as: physical phenomena, chemical phenomena, biological phenomena.

Heat equation will be considered in our study under specific conditions.

Objectives

The aims of this dissertation are:

1. To study the recent literature survey on FTCS method and Adomian Decomposition Method for heat equation.

2. To investigate the performance (in particular accuracy and stability) of these methods when applied to a real problem rather than a simple illustrative problem.

3. To compare between FTCS method and Adomian Decomposition Method under different values.

Overview

  This study consists of six chapters. Chapter one provides the topic to be studied and the aims to be achieved. Chapter two states basis of the finite difference method and some illustrative examples. Chapter three and chapter four we will discuss in greater detail the FTCS method and Adomian decomposition method respectively. Recent literature review on FTCS method and ADM method for heat equation are presented.

  Chapter five presents numerical experiments for one dimensional heat equation using MATLAB software for these methods. Chapter six provides the conclusion and suggestions for further future work.

CHAPTER TWO

Finite Difference Method

Introduction

  Finite difference method (FDM) is widely used for the solution of partial differential equations of heat, mass and momentum transfer. Furthermore, it is very easy to learn and apply to many problems in physics and engineering. Extensive amount of literature exist on the applications of this method for solving such problems. The accuracy of the (FDM) can be easily tested by the order of the truncation error in the Taylor’s series expansion (Ozisik, 1994).

  The idea of the (FDM) proceeds by replacing the derivatives in the partial differential equations with finite difference approximations. This gives large systems of algebraic equations which usually can be solved in the area of partial differential equations. Algebraic equations can readily solved using a computer (LeVeqe, 2007).

  In this chapter we will treat to the Finite difference method and other related topics in more details. Some of what this chapter includes is adapted from the Lecture notes book MAT518 (2008).

Finite Difference Approximations

  Suppose a function u=u(x,y). Divide the first quadrant of the xy plan into uniform rectangles by grid lines parallel to x axis (uniform length) and grid lines parallel to y axis (uniform length ). See Figure (2.1).

A  Assume we have a continues function u, so the Tylor’s series expansion for uxi+x,yj about (xi,yj) is :

uxi+x,yj=uxi,yj+x1!uxxi,yj+x22!uxxxi,yj+x33!uxxxxi,yj+x44!uxxxxxi,yj+…

  Suppose we decide to truncate the series on right hand side (RHD) of (2.1) beginning with the 3rd term. If x is sufficiently small, he 4thA  and higher terms are much smaller than the 3rd term. Then we write:

uxi+x,yj=uxi,yj+x1!uxxi,yj+0x2

 0(x2) Means the sum of the truncated terms (i.e. the truncation error) is in absolute term most a constant multiple ofA  x2 .Divide (2.2) by (x) and after rearranging the terms once we get:

uxxi+x,yj=uxi+x,yj-uxi,yjx+0x 

uxi+x,yj-uxi,yjxA A  Which is called a forward difference approximation to the partial derivative ofA‚u‚x at the points xi,yj and it is known to be the first order accurate or {0(x)} accurate. IfA  x is reduced by fifty-present (50%), then the truncation error is likewise reduced by 50%.

The Taylor’s series expansion for uxi-x,yj about (xi,yj) is :

uxi+x,yj=uxi,yj-x1uxxi,yj+x22uxxxi,yj-x33uxxxxi,yj+x44uxxxxxi,yj+… (2.4)

  Suppose we decide to truncate the series on right hand side (RHS) of equation (2.4) beginning with the 3rd term. If x is sufficiently small, the 4th and higher terms are much smaller than the 3rd term. So we write:

uxi-x,yj=uxi,yj-x1!uxxi,yj+0x2(2.5)

By dividing equation (2.5) by (x) and rearrange the terms once we get:

uxxi,yj=uxi+x,yj-uxi,yjx+0x(2.6)

uxi+x,yj-uxi,yjxA A  Which is called a backward difference approximation to the partial derivative ofA ‚u‚x at the points uxi,yj and it is known to be first order accurate or {0(x)} accurate.

If we subtract equation (2.4) from equation (2.1) we get:

uxi+x,yj- uxi-x,yj=2xuxxi+yj+0x3  (2.7)

By dividing equation (2.7) by (2x) and rearrange the terms once we get:

uxxi,yj=uxi+x,yj-uxi-x,yj2x+0x2  (2.8)

uxi+x,yj-uxi-x,yj2x is called central difference approximation to the partial derivative‚u‚x at the points uxi,yj which is known as second order accurate or {0(x2)} (Lecture note book MAT 518,2008).

  Consider the figure (2.2) below, it is simple that the equation (2.8) approximates the slop of the tangent at point P by the slop of the chord AB, and is said to be central difference approximation. Moreover, we can also approximate the slope of tangent P by the slop of chord PB which gives a forward difference approximation formula. Likewise, the slop of the chord AP gives a backward difference approximation formula (smith, 1985).

  Adding equation (2.1) and equation (2.4) we get a similar approximation for ‚2u‚x2A  at the points uxi,yj

uxi+x,yj+ uxi-x,yj=2uxi,yj+x2uxxxi+yj+0(x4) (2.9)

Dividing (2.9) by (x2)A A A  and rearrange the terms to get :

uxxxi+yj=xi+x,yj-2xi+yj+uxi-x,yjx2+0x2  (2.10)

xi+x,yj-2xi+yj+uxi-x,yjx2A  Which is called the central approximation to the second partial derivative ‚2u‚x2A  and it is the second order {(x2)} accurate.

 Similarly as the previous steps, we can get the approximations to the partial derivates uy and uyy to obtain a forward, backward and central difference approximation y

(Look at Figure 2.3) as follows:

uyxi,yj= uxi,yj+y-uxi,yjyA  (Forward difference approximation)  (2.11)

uyxi,yj=A A  uxi,yj-uxi,yj-yy (Backward difference approximation)  (2.12)

A  uyxi,yj= uxi,yj+y-uxi,yj-y2y(Central difference approximation) (2.13)

uyyxi,yj= uxi,yj+y-2xi,yj+uxi,yj-y2y2A  (Central difference approximation) (2.14)

  Finite difference method is an effective method, but it also suffers from a drawback. A major drawback of the finite difference method appears in its incapability to deal efficiently with the solution of the problems over arbitrarily shaped complex geometries. This happens because of the interpolation difficulties between the interior pints and the boundaries conditions in order to construct the (FDM) expressions for nodes points next the boundaries conditions (Ozisik, 1994).

CTCS Method

  We illustrate the CTCS method (Central Time Central Space) as an example of the finite difference method, this example is taken from (Burden; Faires, 2005).

Example 2.1:

Consider the IBVB below:

2u‚t2-‚2u‚x2=0, A  0<x<1 ,A A  0<tA  ;  (2.15)

With boundary conditions:

u0,t=u1,t=0 ,  0<t;

And initial conditions:

ux,0=sinx  ,  

utx,0=0,  0

  Divide the domain into a set of uniform rectangles with grid lines parallel to the x axis and grid lines parallel to the t axis. Suppose we use x=0.20 (i.e. we use 5 equal sub-intervals), consider the figure2.4 below:

  The grid points are (xi,yi) where i = 1,2,3,4 and xi=ix, and n=1,2,3,4,…A  and tn=nt. The initial conditions satisfies the values of u at points x0,x1,x2,x3,x4,x5 at t = 0. And The boundary conditions satisfies the values of u at points x=x0 and x=x5A  t>0.

We use the central differences forA ‚ 2‚ t2A  and‚ 2u‚x2  once we get:

uin+1-2uin+uin-1t2=ui+1n-2uin+ui-1nx2  (2.16)

uin+1=21-r2uin+r2ui+1n+ui-1n-uin-1  (2.17)

Where r in our example is:

r= A  tx (2.18)

 Equation (2.16) is called CTCS formula, which is an explicit formula.

If values of u at time level(n) and (n+1)

  To implement the stability in CTCS method we make sure that the time step t must be chosen such as which meanstx. We assume that t=0.20 and x=0.20 , so r=1.

Substituting the value of r in the formula (2.17) gives:

uin+1=ui+1n+ui-1n-uin-1  (2.19)

We now study the implementation of CTCS, so we start with:

i=1,n=0

u11=u20+u00-u1-1 (2.20)

  All the values on right hand side are known exceptA  u1-1, so we can compute the value of u1-1A  by using the central difference to the lift hand side of the initial condition utx,0=0, which gives u11=u1-1, substituting this value to equation (2.20)A  yields:

u11= u20+u002  (2.21)

  Substituting the valuesA  u00=0, u20=sin0.4=0.951 to the equation (2.21) we get, u11=0.476, which is the approximation to u(0.2,0.2).

i=2,n=0

u21=u30+u10-u2-1  (2.22)

 Here also all values on right hand side are known exceptA  u2-1, so we can compute the value of u2-1A  by using the central difference to the lift hand side of the initial condition utx,0=0, which gives similarlyA  u21=u2-1, substituting this value to equation (2.22)  gives:

u21=  u30+u102 (2.23)

  Using the values u10=sin0.2=0588, u30=sin0.60.951 in the equation (2.23) we get, u11=0.476, which is the approximation to u(0.4,0.2).

i=3,n=0

u31=u40+u20-u3-1 (2.24)

  All the terms on right hand side are known except  u3-1, so we can compute the value of u2-1  by using the central difference to the lift hand side of the initial condition utx,0=0, which gives similarlyA  u31=u3-1, substituting this value to equation (2.24)  gives:

u11=A  u40+u202 (2.25 Using the values u20=sin0.4=0.951, u40=sin0.8=0.588 in the equation (2.25) we get, u31=0.769, which is the approximation to u(0.6,0.2).

i=4,n=0

u41=u50+u30-u4-1   (2.26)

  All the terms on right hand side are known exceptA  u4-1, so we can compute the value of u4-1A  by using the central difference to the lift hand side of the initial condition utx,0=0, which gives similarlyA  u4-1=u41, substituting this value to equation (2.26)A  gives:

u41=u50+u302  (2.27)

  Using the values u30=sin0.6Ï€=0.951, u50=0, in the equation (2.27)

we get, u41=0.476, which is the approximation to u(0.8,0.2).

By using the pervious values we can easily calculate the second row, so we start with:A 

i=1,n=1

u12=u21+u01-u10=0.769+0-0.88=0.181

i=2,n=1

u22=u31+u11-u20=0.769+0.476=0.951=0.294

i=3,n=1

u32=u41+u21-u30=0.476+0.764-0.951=0.289

i=4,n=1

u42=u51+u31-u40=0+0.769-0.588=0.181

So, the computational continues till final time level required.

Crank –Nicolson method

  Crank –Nicolson method is an example of an implicit method, it can be use for solving diffusion equations (parabolic pde’s), for general, consider the IBVP equation:

‚u‚|i+12=Aµ2‚u2‚x2|ij+12

  We can approximate this equation by using the central difference approximation to ‚u2‚x2 (space derivative) and central difference approximation to ‚u‚tA  (time derivative)

at pointsA  [ix, (j+12)t], in other word, we can approximate the equation at non grid points (xi,tj+1/2) by:

  uij+1-uijt= Aµ2ui+1j+12-2uij+ui-1j+12  (2.28)

RHS of equation (2.28) can be written as:

  Aµ22ui+1j+12-2uij+ui-1j+12x2+ui+1j-2uij+ui-1jx2  (2.29)
Equation (2.28) can be written as:

-rui-1j+1+2-2uij+1-rui+1j+1=rui-1j+2-2uij-rui+1j  (2.30)

  This formula for Crank –Nicolson method and it is convergent and stable for all values of {r} (smith, 1985).

  In Crank –Nicolson method, even we know the values of u at time level n from the IC and BC conditions, uij+1 cannot be calculated directly since ui-1j+1A  andA  ui+1j+1 are known, that is why the Crank –Nicolson method is called an implicit method.

The following figure illustrates Crank –Nicolson method:

The concept of stability

Various of finite difference methods were applied for solving parabolic partial differential equations which forms the diffusion equations. All of these methods essentially reduce the solution of pde’s to algebraic system of equations which can be done using the computer. It also known that the numerical computations carried out, even on a high-tech computer, it can be accurate just for a finite numbers of decimal places, which means, at any stage of calculations, someA  round-offA  error will be introduced. Let E be an error introduced into the calculations at all grid points{i,n,ui-n}. The resulting numerical solution of the differential equation will be:

‚u‚t =β2 ‚2u‚x2A   (2.31)

Let uin be the corresponding solution of the differential equation with finite difference method and suppose there is no round-off error in each calculation. Then we state:

  Ein=uin-ui-n  (2.32)

Stability is very important in getting the numerical solution of partial differential equations by using the finite difference methods. The solution of the finite difference equation is known to be stable, if the error does not increase exponentially as we proceed from one time step to another. In another word, Ein should be bounded as n' for a fixed x and t , or when both x and t tends to zero. There are two main methods to investigate whether a difference method is stable:Von-Neumann Fourier series expansion method and matrix method (Rao,2007).

CHAPTER THREE

FTCS Method

Introduction

  FTCS method (Forward-Time Centered-Space) is one of some important and powerful finite difference methods used for solving parabolic partial differential equations that model many physics phenomenal, such that heat, mass transfer in fluids and unsteady behavior of fluid flow past bodies (Rao,2008).

  For parabolic equations, FTCS approximates the partial derivative ut by the

(first-order Forward-Time) approximation and the partial derivative uxx by the

(second-order Centered-Space) approximation (Hoffman, 2001).

In this chapter, we will treat the FTCS method and other relevant subjects in more general details.

  In our study and specifically in this chapter, FTCS method will be considered and used for solving one dimensional linear heat equations.

  First of all and for more illustrations, let us consider a one-dimensional diffusion equation at the mesh grids points (i,n) as follow (Rao,2008):

‚u‚tA  =α2A  ‚u2‚x2  (3.1)

By introducing Forward difference approximation to the ut and Central difference approximation to the uxx , equation (3.1) becomes:

uxi,tn+t-uxi,tnt=α2uxi+x,tn-2uxi,tn+uxi-x,tnx2 A A  (3.2)

 Notice that, uxi,tn+t,uxi,tn, uxi+x,tnA A  andA A A  uxi-x,tn in equation (3.2) are not actual values of u at  uxi,tn+t,uxi,tn,A  uxi+x,tn

and uxi-x,tn But approximations to u at these points.

So Eq(3.2) can be written as: 

uin+1-uin t=α2ui+1n-2uin+ui-1nx2  (3.3)

So that;

uin+1=uin+α2tx2(ui+1n-2uin+ui-1nA  (3.4)

So we define

  α2tx2=r  (3.5)

Then, equation (3.5) becomes:

  uin+1=uin+r(ui+1n-2uin+ui-1n)ormula, since uin+1 can easily computed if all values of u at time level n are known (Rao, 2008).

Stability analysis for FTCS: Fourier’s method

By formulating the Fourier forms component solutions in the form:

uin=Ï„ne-1A  iθ  (3.7)

Into equation (3.6) and we assume that{A  -1 =I} , yields:

Ï„n+1eIθi=rÏ„neIθ(i-1)+1-2rÏ„neIθi+rÏ„neIθ(i+1)  (3.8)

By dividing Eq (3.8) byA  Ï„neIθi and rearrange the terms we get:

  Ï„=1+reIθ-2+e-Iθ(3.9)

  Ï„2=12-r+r(eIθ-e-Iθ)2  (3.9.1)

Using the identity:

cosθ=(eIθ+e-Iθ)/2,  1-cosθ=2sin2(θ/2)  (3.10)

In Eq (3.9.1)A  and multiply Eq (3.9.1)A A  by 2 gives:

Ï„=14rsin2(θ/2)  (3.10)

For Fourier’s method, a method is stable if and only if Ï„‰¤1 for all values of θ.

FTCSA  method to be stably,A  Fourier method requires:

1-4rsin2(θ/2)‰¤1A A  €Î¸, so that -1‰¤1-4rsin2(θ/2)‰¤1. This inequality is satisfied if r‰¤12 A A  (Kincaid, Cheney, 1996; smith, 1985;Thomas,1999; LeVeqe, 2007).

Stability analysis for FTCS: Matrix method

Concept of matrix

In this field we consider five important concepts of matrix which accedes to prove the the stability requirement for the FTSC by using matrix method:

1st) the norm of the matrix A is a real positive number given a measure of the ‘size’

Of the matrix and must satisfy the following requirements:

1) A>0 if A‰ 0 and A=0 if A=0.

2) cA=cA for a real or complex scalar c.

3) AB‰¤AB.

4) A+B‰¤A+B.

5) A1=maxiai,j.

6) A=maxjai,j.

7) A2=max?¡{λ,λ an eigenvalue of ATA12}.

8) limxâ'Ax=0 if A‰¤1.

9) limx'Ax=0 if and only if λi<1 for all eigenvalue.

2nd) let λi be an eigenvalue of the nA—n matrix A and Xi the corresponding eigenvector, hence

AXi=λiXi.

Therefore λi‰¤A for i=1…,n hence Ï?(A)‰¤A.

3rd) the eigenvalue of a common tridiagonal matrix:

The eigenvalue of nA—n matrix

abA  ‹¯‹¯0c ab‹® c ‹±‹® ‹±‹±‹±‹® ‹±‹±b0‹¯ca

Are λi=a+2bc coskÏ€n+1for k=1,…,n

Where a , b and c may be real or complex number.

4th) Gershgurian’s first theorem:

In this theorem, Gershgurian proved that the largest of the moduli of the eigenvalue of a given square matrix A cannot exed the largest sum of the moduli of the elements along any row or any column, that is to say:

A‰¤A1  or Ï?(A)≤A

5th) Gershgurian’s second theorem:

Let pkbe the sum of the moduli of the elements along kth row excepting the diagonal elements ak,k , then each eigenvalue of matrix A lies inside or on the boundary of the as least one of the circles:

  λ-ak,k=pkA A  (Smith, 1985).

So, we consider the following linear and non-homogeneous heat equation:

ut=uxx+f(x,t)

With the initial condition ux,0=g(x) and boundary conditions as follows:

  u0,t=pt

  ul,t=qt

The explicit FTCS method is used and obtains:

uin+1=sui-1n+1-2suin+sui+1n+sui+1n+kf(ih,nk)A  ,A  i=1,…,M-1A  n=0,1,…,N-1

Thus we have

u1n+1u2n+1‹®‹®‹®uM-1n+1=1-2ss0……0s1-2ss00‹±‹±‹±‹®‹®‹±‹±‹±‹±‹®‹®‹±‹±s0‹¯‹¯‹¯s1-2su1nu2n‹®‹®‹®uM-1n+su0n+kf(h,nk)kf(2h,nk)‹®‹®kf((M-1)h,nk)su0M+kf(Mh,nk)

Thus Un+1=AUn+B.

It is clear that the matrix A can be written as:

1-2ss0……0s1-2ss00‹±‹±‹±‹®‹®‹±‹±‹±‹±‹®‹®‹±‹±s0‹¯‹¯‹¯s1-2s=

10……00010‹®‹±‹®‹®‹±‹®01000……01+s-21……001-2‹±0‹®â±‹±‹±‹®‹®‹±‹±‹±‹®0‹±-2100……1-2

In other words,

A=IM-1+sTM-1

Where matrix IM-1is the identity matrix and matrix TM-1 is a tridiagonal matrix.

Thus, the eigenvalue of TM-1 can be obtained as the following:

λk=-4sin2kÏ€2M, k=0,1,…,M-1.

So the equation will be stable where:

A2=max-4sin2kÏ€2M‰¤1

Thus h'0 , M'  give usA  0<s‰¤12  (smith,1985;Rao, 2007).

Consistency analysis for FTCS: Taylor’s series

To show that the FTCS method is consistent with the Eq (3.1) we substitute the exact solution of the equation ut=α2uxx into the formula (3.6), we abstain:

uxi,tn+1=uxi,tn+ruxi+1,tn-2uxi,tn+uxi-1,tn (3.10)

Expand uxi,tn+1 about xi,tn, we obtain:

uxi,tn+1=uxi,tn+t=

uxi,tn+t‚u‚tx=xi,t=tn+t22‚2u‚t2x=xi,t=tn+t33!‚3u‚t3x=xi,t=tn+t44!‚4u‚t4x=xi,t=tn+…  (3.11)

Expand uxi+1,tn about xi,tn, we obtain:

uxi+1,tn=uxi+x,tn=

uxi,tn+x‚u‚xx=xi,t=tn+x2‚2‚x2x=xi,t=tn+x33‚3u‚x3x=xi,t=tn+x44‚4u‚x4x=xi,t=tn+…  (3.12)

Expand uxi-1,tn about xi,tn, we get:

uxi-1,tn=uxi-x,tn=

uxi,tn-x‚u‚xx=xi,t=tn+x22‚2u‚x2x=xi,t=tn-x33‚3u‚x3x=xi,t=tn+x44‚4u‚x4x=xi,t=tn-…  (3.13)

If the expansions (3.11), (3.12) and (3.13) are substituted into Eq (3.10) and by considering thatA A  r=α2∆t∆x2, we obtain:

uxi,tn+t‚u‚tx=xi,t=tn+t22‚2u‚t2x=xi,t=tn+t33‚3u‚t3x=xi,t=tn+t44‚4u‚t4x=xi,t=tn+…=

uxi,tnα2tx2uxi,tn+x‚u‚xx=xi,t=tn+x2‚2u‚x2x=xi,t=tn+x33‚3u‚x3x=xi,t=tn+x44‚4u‚x4x=xi,t=tn+2uxi,tn+uxi,tn-x‚ˆ‚xx=xi,t=tn+x22‚2u‚x2x=xi,t=tn-x33‚3u‚x3x=xi,t=tn+x44‚4u‚x4x=xi,t=tn-…

Then we get:

‚tx=xi,t=tn+t22‚2u‚t2x=xi,t=tn+t33‚3u‚t3x=xi,t=tn+…

=α2t‚2u‚x2x=xi,t=tn+x212‚4u‚x4x=xi,t=tn+…  (3.14)

By dividing both sides by ∆t and rearrange the terms, then we get:

‚u‚tx=xi,t=tn+t2‚2u‚t2x=xi,t=tn+t23‚3u‚t3x=xi,t=tn+… =

α2‚2‚x2x=xi,t=tn+x212‚4u‚x4x=xi,t=tn+…  (3.15)

As t'0 andA x'0 Eq (3.15)A  approaches ut=α2uxx , which is the diffusion equation, so the FTCS method for Eq(3.15) is consistent with ut=α2uxxA  (Hoffman, 2001).

Convergence

 The convergence of FTCS is related to the Consistency and stability, which is answered by study of both (Consistency & stability).

FTCS is said to be convergence if and only if FTCS is consistent and stable, otherwise it is not convergence (Hoffman, 2001;LeVeqe, 2007).

 

Local truncation error

  By return to the equation (3.15), the FTCS formula for ut=α2uxx is:

  uin+1-uint= α2ui+1n+2uin+ui-1nx2  (3.16)

Let;

Fi,nu=uin+1-uint-α2ui+1n+2uin+ui-1nx2=0  (3.17)

therefore;

Ti,n= Fi,nu=uin+1-uint-α2ui+1n+2uin+ui-1nx2=0  (3.18)

By Taylor’s expansions:

ui+1n=uxi+x,tn= uxi,tn+x‚u‚xx=xi,t=tn+x2!‚2u‚x2x=xi,t=tn+x33‚3u‚x3x=xi,t=tn+x44‚4u‚x4x=xi,t=tn+…(3.19)

ui-1n=uxi-x,tn=uxi,tn-x‚u‚xx=xi,t=tn+x22‚2u‚x2x=xi,t=tn-x33!‚3u‚x3x=xi,t=tn+x44!‚4u‚x4x=xi,t=tn-… (3.20)

uin+1=uxi,tn+t=uxi,tn+t‚u‚tx=xi,t=tn+t2!‚2‚t2x=xi,t=tn+t33!‚3u‚t3x=xi,t=tn+t44!‚4u‚t4x=xi,t=tn+… (3.21)

Since ;

  ‚u‚t-α2‚u2‚x2= 0 

Substituting Eq(s): (3.19),(3.20) and(3.21) in Eq(3.18), once we get:

Ti,n=‚u‚u-α2‚u2‚x2x=xi,t=tn+t2‚2u‚t2x=xi,t=tn-α2x212‚4u‚x4x=xi,t=tn+t26‚3u‚t3x=xi,t=tn-α2x4360‚6u‚x6x=xi,t=tn+…  (3.22)

Therefore, the principal part of the local truncation error is:

t2‚2u‚t2x=xi,t=tn-α2x212‚4u‚x4x=xi,t=tn+…  (3.23)

A A  So, this is the local truncation error for the FTCS method for the diffusion equation at the points(ix,nt). And local truncation error for FTCS method for the diffusion equation is said to be 0t+0(x)2 (smith, 1985).

Problem statement

First problem:

  Consider the second order one dimensional linear heat equation:

utx,t=Ï(x,t)uxxx,t=x2tuxxx,t  on 0,1A—0,T  (3.24)

Subjected to the boundary conditions:

u0,t=0.

u1,t=exp(t2).

And the initial condition:

u(x,0)=fx=x2

The exact solution for this problem is given by (A.Soufyan; M.Boulmalf, 2005):

ux,t=x2exp(t2)

FTCS method for first problem:

  By applying the Forward difference approximation to ut(x,t) andA  Central difference approximation toA  uxx(x,t), so Eq(3.24) becomes:

uxi,tn+t-u(xi,tn†t=ix2(nt)uxi+x,tn-2uxi,tn+uxi-x,tnx2

Multiplying both sides by t(x2) and rearrange the terms we obtain:

x2uin+1-uin=ix2nt(t)ui+1n-2uin+ui-1n

So:

x2uin+1=x2uin+ix2nt(t)ui+1n-2uin+ui-1n

Dividing both sides by (x2) yields:

uin+1=uin+tx2ix2nt(t)ui+1n-2uin+ui-1n

Let r=tx2

So, the FTCS formula for equation 3.24 is:

uin+1=uin+rix2nttui+1n-2uin+ui-1n (3.25)

Second problem:

  Consider the second order one dimensional linear heat equation:

utx,t=uxxx,t+Ïx,t=uxxx,t+exp-x(cos t-sint)  (3.26)

Subjected to the boundary conditions:

u0,t=sint

u1,t=sin(t/e)+1

And the initial condition:

ux,0=fx=x

The exact solution for this problem is given by (Bongsoo,2007):

ux,t=x+exp(-x)sint)

FTCS method for second problem:

  By applying the Forward difference approximation to ut(x,t) andA  Central difference approximation toA  uxx(x,t), ),, so Eq(3.26) becomes:

uxi,tn+t-u(xi,tn†t=uxi+x,tn-2uxi,tn+uxi-x,tnx2+e-ix(cosnt)-sinnt)

Multiplying both sides by t(x2) and rearrange the terms and by assuming that

μ=e-ix(cosnt)-sinnt) so we get:

x2uin+1-uin=tui+1n-2uin+ui-1n+(tx2)μ

So it becomes:

x2uin+1=x2uin+tui+1n-2uin+ui-1n+(tx2)μ

By solving the equivalent expression by uin+1 we get:
uin+1=uin+tx2ui+1n-2tx2uin+tx2ui-1n+(t)μ

Let r=tx2A  so:

uin+1=uin+rui+1n-2ruin+rui-1n+(t)μ

Substituting the value of μ and rearrange the terms once we get:

uin+1=rui-1n+1-2ruin+rui+1n+te-ix(cosnt-sinnt)  (3.27)

Which is the FTCS formula for equation 3.26.

Third problem:

Consider the second order one dimensional linear heat equation:

utx,t=uxxx,t+Ï•x,t=uxxx,t+exp-x(cosh t-sinht)  (3.28)

Subjected to the boundary conditions:

u0,t=sinht

u1,t=e.sinht+t+1/6

And the initial condition:

ux,0=fx=x3/6

The exact solution for this problem is given by (Bongsoo,2007):

ux,t=exp(x)sinht+x3/6+xt

FTCS method for third problem:

  By applying the Forward difference approximation to ut(x,t) andA  Central difference approximation toA  uxx(x,t), so Eq(3.27) becomes:

uxi,tn+t-u(xi,tn)t=uxi+x,tn-2uxi,tn+uxi-x,tnx2+e-ix(coshnt)-sinhnt)

Let β=e-ix(coshnt)-sinhnt)

Multiplying both sides by t(x2) and rearrange the terms and by assuming that

μ=e-ix(cosnt)-sinnt) so we get:

x2uin+1-uin=tui+1n-2uin+ui-1n+(tx2)β

So it becomes:

x2uin+1=x2uin+tui+1n-2uin+ui-1n+(tx2)β

By solving the equivalent expression by uin+1 we get:
uin+1=uin+tx2ui+1n-2tx2uin+tx2ui-1n+(t)β

Let r=tx2A  so:

uin+1=uin+rui+1n-2ruin+rui-1n+(t)β

Substituting the value of β and rearrange the terms ,the FTCS formula for Eq(3.28):

uin+1=rui-1n+1-2ruin+rui+1n+te-ix(coshnt-sinhnt)  (3.29)

CHAPTER FOUR

ADM Method

Introduction

  Over the last 10 years, many numerical methods are aimed to solve non-linear partial differential equations have appeared in many researcher’s studies. However, most of these methods require a tedious analysis work which makes unnecessary linearization, perturbations, restrictive methods and assumptions which may change the problem itself being solved sometimes. Adomian decomposition method (ADM) overcame the major part of these difficulties were suffered and it was successfully applied in many applications in applied sciences (JAVIDI, M. G., A, 2007; LUO, X.-G, 2005).

  In our study, Adomian decomposition method (ADM) will be considered and used for solving one dimensional linear heat equation.

Literature survey

  In the early 1980s, George Adomian developed a fruitful powerful method for the computational of exact solution for linear and non-linear partial differential equations (pde’s) (Adomian, G, 1984; Adomian, G, 1988; Adomian, G, 1990; Adomian, G, 1998).

 wide class of problems in mathematics, physics, engineering, biology, chemistry and in other sciences are presented and solved using Adomian decomposition method by many authors (Adomian, G, 1984; Adomian, G, 1998; Wazwaz, A.-M. E.-S., Salah, M, 2001; Wazwaz, A.-M, 2000; Chen, W. L., Zhengyi, 2004; Jang, B,2007; Amani, R., A.Sadeghi,J ,2007).

  Furthermore, the method provides an effective tools for analytical solutions of a wide class of dynamical systems representing real physical and engineering problems, which generally converge very rapidly in order to reach an accurate solution, is needs less calculations to solve parabolic partial differential equations comparison with other methods (Adomian, G,1984; Adomian, G,1990; Soufyane, A. B., M,2005; Ray, S., S.Bera,K,R,2005; Bildik, N. B., Hatice,2005).

Analysis of the method

  In this section and for convenience of the readers, we review the basic principles of the standard Adomian decomposition method by the following procedures:

  We consider the following differential equation (G.Adomian, 1988):

Lu+Ru+Nu=g(x,t)  (4.1)

Where L represents the highest order derivative which assumed to be invertible, R is the linear differential operator of order less than L, where g(x,t) is the source term, and finally, Nu represents the non-linear terms.

Solving Lu from Eq(4.1), we obtain:

Lu=gx,t-Ru-Nu    (4.2)

Since L is :

L-1Lu=L-1(gx,t)-L-1(Ru)-L-1(Nu)  (4.3)

If  L-1 is second order operator, so L-1 is a two-fold integration operator assumed to be stated from (to=0)  to  t  .

Therefore, Lt , Lx are differential operators defined as:

Lt=‚2‚t2  , Lx=‚2‚x2  and  L-1=0t0t.dtdt .

Solving Eq(4.3) and using the given conditions we obtain:

u=ω+L-1(gx,t)-L-1(Ru)-L-1(Nu)  (4.4)

  Where the function (ω=A+Bt) represents all the terms arising from the given conditions, all the conditions are assumed to be prescribed for a problem aimed to be solved.

  For non-linear equations, the non-linear term Nu=F(u) will be decomposed by an infinite series of the so-called Adomian polynomials, so:

Fu=n=0An

WhereA A A A  An=A  1ndndλnNn=0unλnλ=0 [A  5]

Where An are generated for all kinds of nonlinearity, so Ao depends on u0 only, A1 depends on u0 and u1 ,and so on until An depends on u0 , u1 , ……,un .

For more details about formulating Adomian polynomials, specific algorithms were set in .

The standard Adomian decomposition method decomposes the solution u(x,t) into an infinite series defined by:

ux,t=n=0un(x,t)  (4.5)

Where u0 , u1 , ……,un are usually determined as:

u0=ψ+L-1(gx,t)

u1=-L-1Ru0-L-1(A0)

un+1=-L-1Run-L-1(An) n0  (4.6)

Substituting the values of An in Eq(4.6) leads to the determination of the solution of u in series form defined by Eq(4.5) immediately.

So, the approximate solution using n-terms approximation to u(x) is given by φn such that:

ux,t=limn'Ïn

Where Ïn is defined as:

Ïn =k=0n-1uk(x,t)

Which can be written as:

Ï1 =u0,

Ï2 =u0+u1,

Ï3 =u0+u1+u2,

Ï4 =u0+u1+u2+u3,

Ïn =u0+u1+u2+u3+…+nn-1.

So , the absolute error is defined as:

absolute error =exact solution-approximate solution

  Abdul-Majid Wazwaz have investigate the existence of the so-called noise terms ’self–cancelling’ phenomenon, which if it exists, it gives a useful tool to fast convergence of the solution by the first two iterations only (u0 & u1). It is noted that the noise terms phenomenon may appear only in non-homogenous problems, whereas it is not exist in homogenous problems. The remaining “non-canceling”? terms may give the exact solution for the non-homogenous problem .

  In our study, the linear operators of the Adomian decomposition method (ADM) will be considered and used for solving one dimensional linear heat equations.

Problem statement

First problem:

Consider the first problem which presented in chapter three and we recall it as:

utx,t=Ï(x,t)uxxx,t=x2tuxxx,t   on 0,1A—0,  (4.7)

Subjected to the boundary conditions:

u0,t=0.

u1,t=exp(t2).

And the initial condition:

u(x,0)=x2

The exact solution for this problem is given by (A.Soufyan; M.Boulmalf, 2005):

ux,t=x2exp(t2)

ADM for first problem:

By using the standard Adomian decomposing method to determine the individual terms from Eq(4.7), so Eq(4.7) becomes:

Ltux,t=x2tLxux,t  4.8

Since Lt is investable, so we have:

Lt-1Ltux,t=Lt-1(x2tLxux,t) 4.9

By using the given initial condition, thus we have:

ux,t=ux,0+Lt-1(x2tLxux,t)  4.10

Substituting the series 4.5 into (4.10) the we obtain:

u0=x2

u1=x20tt‚2u0‚x2dt=x2t2,

u2=x20tt‚2u1‚x2dt=12!x2t4,

u3=x20tt‚2u2‚x2dt=13!x2t6,

u4=x20tt‚2u3‚x2dt=14x2t8,

u5=x20tt‚2u‚x2dt=15x2t10,

u6=x20tt‚2u5‚x2dt=16x2t12,

u7=x20tt‚2u6‚x2dt=17x2t14,

un=x20tt‚2un-1‚x2dt=1n!x2t2n,

And so on, the exact solution for the given problem is in series form given by:

ux,t=x2exp(t2)

Which easily can be verified.

Second problem:

by Considering the second problem which presented in chapter three and we recall it as:

utx,t=uxxx,t+φx,t=uxxx,t+exp-x(cos t-sint) (4.11)

Subjected to the boundary conditions:

u0,t=sint

u1,t=sin(t/e)+1

And the initial condition:

ux,0=fx=x

The exact solution for this problem is given by (Bongsoo,2007):

ux,t=x+ex¡(-x)sint)

ADM for second problem:

Applying the ADM for this problem yields:

Ltux,t=Lxux,t+exp-xcos t-sint  (4.12)

Reapplying the same procedures in the previous problem then we obtain:

ux,t=ux,0+Lt-1(Lxux,t)+Lt-1(exp-xcos t-sint)  (4.13)

So, the individual terms are:

u0=x+exp-x(sint+cost-1),

u1=exp-x(1-cost+sint-t),

u2=exp-x(1+t-12t2-sint-cost),

u3=exp-x(-1+t+12t2-13!t3+cost-sint),

u4=exp-x(-1-t+12!t2+13!t3-14!t4+sint+cost),

u5=exp-x(1-t-12!t2+13!t3+14!t4-15!t5-cost+sint),

u6=exp-x(1+t-12!t2-13!t3+15!t5-16!t6-sint-cost),

u7=exp-x(-1+t+12!t2-13!t3-14!t4+16!t6-17!t7+cost-sint),

u8=exp-x(-1-t+12!t2+13!t3-14!t4-15!t5+17!t7-18!t8+sint+cost),

Third problem:

Finally, by considering the third problem which presented in chapter three and we recall it as:

utx,t=uxxx,t+Ï•x,t=uxxx,t+exp-x(cosh t-sinht)  (4.14)

Subjected to the boundary conditions:

u0,t=sinht

u1,t=e.sinht+t+1/6

And the initial condition:

ux,0=fx=x3/6

The exact solution for this problem is given by (Bongsoo,2007):

ux,t=exp(x)sinht+x3/6+xt

ADM for third problem:

Applying the ADM for this problem yields:

Ltux,t=Lxux,t+exp-xcosh t-sinht  (4.15)

Reapplying the same procedures in the previous problems then we obtain:

ux,t=ux,0+Lt-1Lxux,t+Lt-1(exp-xcosh t-sinht)  (4.16)

So, the individual terms are:

u0=16x3+expx(sinht-cosht+1),

u1=xt+expx(cosht-sinht-1+t),

u2=expx(sinht-cosht+1-t+12!t2),

u3=expx(cosht-sinht-1+t-12!t2+13!t3),

u4=expx(sinht-cosht+1-t+12!t2-13!t3+14!t4),

u5=expx(cosht-sinht-1+t-12!t2+13!t3-14!t4+15!t5),

u6=expx(sinht-cosht+1-t+12!t2-13!t3+14!t4-15!t5+16!t6),

u7=expx(cosht-sinht-1+t-12!t2+13!t3-14!t4+15!t5-16!t6+17!t7),

u8=expx(sinht-cosht+1-t+12!t2-13!t3+14!t4-15!t5+16!t6-17!t7+18!t8),

Cite this page