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:

Fx,y,u,∂u∂x,∂u∂y,∂2u∂x2,∂2u∂y2,∂u∂x∂y,…=0

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:  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:

∂u∂t=α∂2u∂x2         ,       Â  for        0≤x≤xf  ,  0≤t≤T

   Â  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 ∆y) and grid lines parallel to y axis (uniform length ∆x). See Figure (2.1).

  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+…                                    (2.1)

  Â  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 4th  and higher terms are much smaller than the 3rd term. Then we write:

uxi+∆x,yj=uxi,yj+∆x1!uxxi,yj+0∆x2                        (2.2)

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

uxxi+∆x,yj=uxi+∆x,yj-uxi,yj∆x+0∆x                     Â  (2.3)

uxi+∆x,yj-uxi,yj∆x  Which is called a forward difference approximation to the partial derivative of  ∂u∂x at the points xi,yj and it is known to be the first order accurate or {0(∆x)} accurate. If  ∆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-∆x1!uxxi,yj+∆x22!uxxxi,yj-∆x33!uxxxxi,yj+∆x44!uxxxxxi,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+0∆x2                        (2.5)

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

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

uxi+∆x,yj-uxi,yj∆x  Which is called a backward difference approximation to the partial derivative of  ∂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=2∆xuxxi+yj+0∆x3              Â  (2.7)

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

uxxi,yj=uxi+∆x,yj-uxi-∆x,yj2∆x+0∆x2                    Â  (2.8)

uxi+∆x,yj-uxi-∆x,yj2∆x 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∂x2  at the points uxi,yj

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

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

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

xi+∆x,yj-2xi+yj+uxi-∆x,yj∆x2  Which is called the central approximation to the second partial derivative ∂2u∂x2  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,yj∆y  (Forward difference approximation)   Â  (2.11)

uyxi,yj=  uxi,yj-uxi,yj-∆y∆y (Backward difference approximation)   (2.12)

  uyxi,yj= uxi,yj+∆y-uxi,yj-∆y2∆y(Central difference approximation)  (2.13)

uyyxi,yj= uxi,yj+∆y-2xi,yj+uxi,yj-∆y2∆y2  (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,   0

With boundary conditions:

u0,t=u1,t=0  ,  0

And initial conditions:

ux,0=sinÏ€x  ,  0≤x≤1 ,

utx,0=0,  0≤x≤1

    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=i∆x, and n=1,2,3,4,…  and tn=n∆t. 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=x5  ∀t>0.

We use the central differences for  ∂2u∂t2  and ∂2u∂x2  once we get:

uin+1-2uin+uin-1∆t2=ui+1n-2uin+ui-1n∆x2                               Â  (2.16)

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

Where r in our example is:

r=   ∆t∆x                                                     Â      (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 r≤1 which means∆t≤∆x. 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 except  u1-1, so we can compute the value of u1-1  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)  yields:

u11= u20+u002                                                    Â  (2.21)

  Â  Substituting the values  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 except  u2-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 similarly  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.6Ï€=0.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 similarly  u31=u3-1, substituting this value to equation (2.24)  gives:

u11=  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 except  u4-1, so we can compute the value of u4-1  by using the central difference to the lift hand side of the initial condition utx,0=0, which gives similarly  u4-1=u41, substituting this value to equation (2.26)  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: 

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∂t|ij+12=µ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∂t  (time derivative)

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

     Â  uij+1-uij∆t= µ2ui+1j+12-2uij+ui-1j+12∆x2                                  (2.28)

RHS of equation (2.28) can be written as:

           Â  µ22ui+1j+12-2uij+ui-1j+12∆x2+ui+1j-2uij+ui-1j∆x2                           Â  (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+1  and  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, some  round-off  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∂x2             Â                 (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∂t  =α2  ∂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,tn∆t=α2uxi+∆x,tn-2uxi,tn+uxi-∆x,tn∆x2   (3.2)

    Notice that, uxi,tn+∆t,uxi,tn, uxi+∆x,tn  and Â  uxi-∆x,tn in equation (3.2) are not actual values of u at  Â  uxi,tn+∆t,uxi,tn,  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-1n∆x2                                     (3.3)

So that;

uin+1=uin+α2∆t∆x2(ui+1n-2uin+ui-1n                              Â  (3.4)

So we define

     Â         α2∆t∆x2=r                   Â              Â      Â  (3.5)

Then, equation (3.5) becomes:

   uin+1=uin+r(ui+1n-2uin+ui-1n)                       Â    Â  (3.6)

Eq(3.6) is the FTCS formula, 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-1  iθ                                             Â       Â  (3.7)

Into equation (3.6) and we assume that{  -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) by  τ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)  and multiply Eq (3.9.1)  by 2 gives:

Ï„=1-4rsin2(θ/2)                                    Â          Â  (3.10)

For Fourier’s method, a method is stable if and only if τ≤1 for all values of θ.

FTCS  method to be stably,  Fourier method requires:

1-4rsin2(θ/2)≤1  ∀θ, so that -1≤1-4rsin2(θ/2)≤1. This inequality is satisfied if r≤12   (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 n×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 n×n matrix

a    Â  b  ⋯⋯0  c    Â  a    b  Â  â‹®     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=pk  (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)  ,  i=1,…,M-1  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 us  0

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!∂2u∂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 that  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α2∆t∆x2uxi,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∂u∂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:

∆t∂u∂tx=xi,t=tn+∆t22!∂2u∂t2x=xi,t=tn+∆t33!∂3u∂t3x=xi,t=tn+…

=α2∆t∂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∂2u∂x2x=xi,t=tn+∆x212∂4u∂x4x=xi,t=tn+…                         Â  (3.15)

As ∆t→0 and  ∆x→0 Eq (3.15)  approaches ut=α2uxx , which is the diffusion equation, so the FTCS method for Eq(3.15) is consistent with ut=α2uxx  (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-uin∆t= α2ui+1n+2uin+ui-1n∆x2                         Â  (3.16)

Let;

Fi,nu=uin+1-uin∆t-α2ui+1n+2uin+ui-1n∆x2=0                    (3.17)

therefore;

Ti,n= Fi,nu=uin+1-uin∆t-α2ui+1n+2uin+ui-1n∆x2=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!∂2u∂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-α2∆x212∂4u∂x4x=xi,t=tn+∆t26∂3u∂t3x=xi,t=tn-α2∆x4360∂6u∂x6x=xi,t=tn+…         (3.22)

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

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

  So, this is the local truncation error for the FTCS method for the diffusion equation at the points(i∆x,n∆t). And local truncation error for FTCS method for the diffusion equation is said to be 0∆t+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,1×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) and  Central difference approximation to  uxx(x,t), so Eq(3.24) becomes:

uxi,tn+∆t-u(xi,tn)∆t=i∆x2(n∆t)uxi+∆x,tn-2uxi,tn+uxi-∆x,tn∆x2

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

∆x2uin+1-uin=i∆x2n∆t(∆t)ui+1n-2uin+ui-1n

So:

∆x2uin+1=∆x2uin+i∆x2n∆t(∆t)ui+1n-2uin+ui-1n

Dividing both sides by (∆x2) yields:

uin+1=uin+∆t∆x2i∆x2n∆t(∆t)ui+1n-2uin+ui-1n

Let r=∆t∆x2

So, the FTCS formula for equation 3.24 is:

uin+1=uin+ri∆x2n∆t∆tui+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) and  Central difference approximation to  uxx(x,t), ),, so Eq(3.26) becomes:

uxi,tn+∆t-u(xi,tn)∆t=uxi+∆x,tn-2uxi,tn+uxi-∆x,tn∆x2+e-i∆x(cosn∆t)-sinn∆t)

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

μ=e-i∆x(cosn∆t)-sinn∆t) so we get:

∆x2uin+1-uin=∆tui+1n-2uin+ui-1n+(∆t∆x2)μ

So it becomes:

∆x2uin+1=∆x2uin+∆tui+1n-2uin+ui-1n+(∆t∆x2)μ

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

Let r=∆t∆x2  so: