"5 - 1Chapter – 5 Numerical Differentiation5.1Introduction to Finite DifferencesNumerical differentiation is method of obtaining the derivation of a function using a numerical technique.In numerical analysis, we get the result in numerical form by computing methods of given data. The baseof numerical analysis is calculus of finite difference which deals with the changes in the dependentvariable due to changes in the independent variable.The finite difference method is used to solve ordinary differential equations that have conditions imposedon the boundary rather than at the initial point.These problems are called boundary-value problems.Suppose the function y = f(x) has the values ? , ? , ? ,……. ? for the equally spaced values? , ? + h, 0 1 2 ? 0 0 ? +2h, …….? + ? h. If y = f (x) be any function then the value of the independent variable ‘x’ is called2 0 argument and corresponding value of dependent variable y is called entry. To determine the value of y?? and for some intermediate values of x is based on the principle of finite difference. This requires three?? types of differences.? Forward Differences ? Central Differences? Backward DifferencesForward Differences The differences? ? , ? ?? ? are called the first forward differences of the function1 - 0 2 - 1,………. ? - ? -1 y = f (x) and we denote these difference by ?? ,? ? .....?? , respectively, where ? is called the0 1 ? descending or forward difference operator.In general, the first forward differences is defined by?? =? -?? ? +1 ? The differences of the first forward differences are called the second forward differences and denoted by2 2 ? ? , ? y etc. Therefore, we have 0 2 ? ? = ? [ ? - ? ] 0 1 0 = ?? - ??10= ( ? - ? ) - ( ? - ? ) 21 1 0=? - 2 ? +?21 0 2 Similarly ? ? = ? [ ? - ? ] 1 2 1 = ?? - ??21 = ( ? - ? ) - ( ? - ? ) 32 2 1=? - 2 ? +?32 1 2 In general we have ? ? = ?? - ??x x+1 xAgain, the differences of second forward differences are called third forward differences and denoted by3 3 ? ? , ? ? etc.0 1 thIn general, the n forward difference is given by n n-1 n-1 ? ? = ?? - ??x x+1 xForward Difference Table 5 - 2Argument Entry FirstSecondThirdForthDifferenceDifferences Differences Differencesy =2 3 4 x ?? ? ? ? ? ? ?f(x)??0 0? - ? = ??1 0 0 2 ? + h ?? ? - ?? = ? ?0 1 10 0 2 2? - ? = ? ?2 1 1 ? ? - ? ? =1 03 ? ? 0 3 3? ? - ? 1 2 4? ? - ? ? = ? ? ? + 2h ? ? = ? ?0 2 21 1 0 0 2 2? - ? = ? ?3 2 2? ? - ? ? =2 13 ? ? 1 2 ? + 3h ?? ? - ?? = ? ?0 3 32 1? - ? = ??4 3 3? + 4h ?0 4 Central Differences The differences? ? = ? ? 1/2, ? ? = ? ? 3/2,? ? = ? ? n-1/2, are called central differences1 - 0 2 - 1 ? - ? -1and d is called central difference operator.2 Similarly-d y 3/2 1/2 = 1 2- = d y and5/2 3/2 22 2 3d y - d y = d y and so on.2 1 3/2Central Difference TableFourth Argument Entry First Differences Second ThirdDifferencesDifferences Differences 2 3 4 x y = f(x) ? ? ? ? ? ? d ???0 0d y1/2 2d y? ? 1 1 1 3d y dy3/2 3/2 2 4 ? ?d ydy2 2 2 2 3d y dy5/2 5/2 2d y ? ? 33 3 d y 7/2? ?4 4 ?? ?? ?? ??5 - 3The working expression of central difference approximation to the first four difference of y is:i y ’ = 1/2h(y – y )i i-1 i+1 2 y ’’ = 1/h (y – 2y + y )i i+1 i i-1 3 y ’’’ = 1/2h (y – 2y + 2y – y )i i+2 i+1 i-1 i-2 2 y ’’’’ = 1/h ((y – 4y + 6y - 4y + y )i i+2 i+1 i i-1 i-2 Backward Differences The differences? ? , ? ?? ? are called the first backward differences of the function1 - 0 2 - 1,………. ? - ? -1 y = f (x) and we denote these difference by?? , ? ? .....? ? , respectively, where ? is called the0 1 ? ascending or backward difference operator.In general, the first backward differences is defined by ?? =? -?? ? ? -1 The differences of the first backward differences are called the second backward differences and denoted 2 2 by ? ? , ? y etc. Therefore, we have 3 22 ? ? = ? [ ? - ? ] 2 2 1 = ?? - ??2 1= ( ? - ? ) - ( ? - ? ) 21 1 0=? - 2 ? +?21 02 In general we have ? ? = ?? - ??x x x-1Again, the differences of second backward differences are called third backward differences and denoted 3 3 by ? ? , ? ? etc.3 4 thIn general, the n backward difference is given bynn-1 n-1 ? ? = ?? - ??x x x-1 Backward Difference TableFourthFirst Second ThirdDifferencesArgument EntryDifferences Differences Differences 2 3 4 x y = f(x) ? ? ? ? ? ? ? ???0 0? – ? = 1 0 ? ?1? ? - ? ? =21 ? + h ?0 1 2? ?22 2? - ? = ? ? - ? ? 3 2 21 3 ? ? = ? ?3 2 3 3? ? - ? ? =? ? - ? ? = 3 2 4 3 4 ? + 2h ?? ? 2 0 2 4 ? ?32 2? - ? = ? ? - ? ? 3 2 4 3 3 = ? ?? ?3 4? ? - ? ? = 43 ? + 3h ?0 32 ? ?4 ? - ? = 43 5 - 4 ? ?4 ? + 4h ?0 4 5.1.1Solved ExamplesExample 1: Form a forward difference table using below data.X: 10 11 12 13 14 15 16Y: 7.989 8.403 8.781 9.129 9.451 9.750 10.031Solution:First, we form forward difference table2 3 4 5 6 x y = f(x) ?? ? ? ? ? ? ? ? ? ? ?10 7.9890.414 11 8.403- 0.0360.378- 0.002 12 8.781- 0.0380.030.3400.028- 0.079 13 9.121- 0.010- 0.0490.1620.330- 0.021 0.083 14 9.451- 0.031 0.0340.299 0.013 15 9.750- 0.0180.281 16 10.0314 Example 2: If U0 = 3, U1 = 12, U2 = 81, U3 = 2000, U4 = 100, Calculate ? U0.Solution:First, we form forward difference table2 3 4 x U = f(x) ?U ? U ? U ? U0 39 1 1260691790 2 811850- 74591919- 5669 3 2000 - 3859-1900 4 100 5 - 54 Hence ? U = - 74590 Example 3: Construct the table of difference for the data below:x:0 1 2 3 4 y: 1.0 1.5 2.2 3.1 4.64 Evaluate ? f(2).Solution:First, we form forward difference table2 3 4 x y= f(x) ?? ? ? ? ? ? ?0 1.00.5 1 1.50.20.70 2 2.20.20.40.90.4 3 3.10.61.5 4 4.64Hence ? f(2) = 0.4.3 Example 4: Show that ? y = y - 3y + 3y - y i i+3 i+2 i+1 i Solution:Step 1 Find ?yi = ?yy -yi i+1 i 2 Step 2 Find ? y 2 i ? y =? ( ? y )i i = ? (y – y )i+1 i= ?y – ?yi+1i = (y – y ) – (y – y )i+2 i+1 i+1 i= yi+2 –2 yi+1+ yi 3 Step 3 Find ? y3 2 i ? y =? ( ? ? )i ? = ? ( y –2 y + y )i+2 i+1i= ? y –2 ? y +? yi+2 i+1i = (y – y ) –2 ( y – y ) + ( y – y )i+3 i+2 i+2 i+1 i+1 i =y – y –2y +2 y +y – yi+3 i+2 i+2i+1 i+1 i=y –3 y +3 y – yi+3 i+2 i+1 i3 Hence Proved ? y = y - 3y + 3y - y i i+3 i+2 i+1 i th th Example 5: What will be the values of the (n+1) differences for a polynomial of n degree?Solution: th 2Consider the n degree polynomial f(x) = A + A x + A x +..... + A x where A , A , A , ..... A are0 1 2 n 0 1 2 nconstants and n is a positive integer.By the definition, we have? f(x) =f(x + h) – f(x)2 n= [A + A (x + h) + A (x+ h) +….+ A (x+ h) ] –0 1 2 n5 - 62 n [A + A x + A x+ ..... + A x ]0 1 2 n2 2 n n = A h + A [(x+ h) – x ]+….+ A [(x+ h) - x ]1 2n2 2 2 2 n n n-1n = A h + A [x + C x h + h - x ]+….+ A [x + C x h + C 1 21 n1 2n-2 2 n n n x h +……..+ C h - x ]n2 n-2 n-1=B + B x + B x +……+B x + n A hx-------------------(1)1 2 3 n-1 n Where B ,B ,B are constants1 2 n-1From equation (1), we see that the first difference of a polynomial of degree n is again a polynomial ofdegree (n – 1).2 ? f(x) =?f(x + h) – ?f(x)2 n -1 = [B + B (x + h) + B (x+ h) +….+n A h (x+ h) ] –1 2 3 n 2 n-2 n-1 [B + B x + B x+ ..... + B x +n A h x ]1 2 3 n-1 n2 2 n -2 n-2 = B h +B [(x+ h) – x ]+….+B [(x+ h) - x ] + 2 3 n-1 n -1 n-1 n A h[(x+ h) - x ]n 2 2 2 2 n-2 n-2 n-3 = B h +B [x + C x h + h - x ]+….+B [x + C x h + 2 3 1 n-11n-2 n-3 2 n-2 n-2 n-2 n -1 n-1 n-2 C x h +……..+ C h - x ] + nA h [x + C x h + 2 n-2 n1n-1 n-3 2 n-1 n n-1 C x h +……..+ C h - x ]2 n-12 n-3 2 n-2 =C + C x + C x +……+C x + n(n-1)h A hx----------(2)2 3 4 n-1 n WhereC2 , C3 ……. ,Cn-1 are constantsFrom equation (2) we see that the second difference of a polynomial of degree n is again a polynomial ofthdegree (n –2). Proceeding in the same way, we will get a zero degree polynomial for the n difference i.e. n n n – n ? f(x) = n (n – 1) (n – 2).... 1 h a xn n= n ! h an th Thus, n difference is constantn+1 nNow ? f(x) = ? [ ? f(x)] n= ? [n! ha ] n= 0 [Since ? C = 0]th th Hence the values of the (n+1) differences for a polynomial of n degree will be zero.2 2 Example 6: Evaluate ? ((5x + 12) / (x + 5x + 6))Solution:(5x + 12)=2 2 3? 2 ? { +}2 (x + 5x + 6) (? + 2) (? + 3) =2 3?{? [] +? []}( ) ( ) ? + 2 ? + 3 =1 1 1 1?{2 [- ] +3 [- ]}(? + 3) (? + 2) (? + 4) (? + 3) =1 1-2? [] - 3? []( ) ( ) (? +2) ? + 3 ? + 3 (? +4) =1 1-2[ - ]-( ) ( ) ? + 3 (? +4) (? +2) ? + 3 1 13[ -]( ) ( ) ? + 4 (? +5) ? + 3 (? +4) 4 6 = +(? +2)(? + 3)(? +4) (? + 3)(? +4)(? +5)5 - 7= 2(5? +16) ( )( )( ) ? +2 ? + 3 ? +4 (? +5) Example 7: Using second order central difference to compute maximum deflection in a simply supportedbeam of span ‘L’ caring a uniform distributed load ‘w’ k N/m. The textural rigidity of the beam isconstant throughout. Use four sub intervals.Solution:?? = ????? / EI?? ? 2? =2 ? 3? = ? h??? ?????3 ? 4? = -??? 4 w per unit lengthi = 1i = 0 i = 2 i = 3 i = 4x = L/4x = 0 1 x = L/2 x = 3L/4 x = L0 2 3 4 y 1y = 0 yyy = 00 2 3 4 For i = 1, h = L/416EI wLx1 wx12 [ ] y2 - 2y1 + y0 = -L2 2 2 16EI w * L * L w * L * L [y2 - 2y1] = *L2 2 * 4 2 * 16 * EI ???? [ ] - =? * ????? For i = 2, h = L/216EI wLx2 wx22 [ ] y3 - 2y2 + y1 = -L2 2 2 16EI w * L * L w * L * L [y3 - 2y2 + y1] = -L2 4 2 * 4 wL4 [ ] y3 - 2y2 + y1 =128EI ??? ?? ?? ?? ?? ??5 - 8For i = 3, h = 3L/416EI wLx3 wx32 [y4 - 2y3 + y2] = -L2 2 2 3wL4 [ ] -2y3 + y2 =512EI wL4 =128EI 5.2 Numerical Solution of Partial Differential EquationA number of physical processes can be described by a partial differential equation (PDE) or a system ofpartial differential equations.For example, conduction of heat in a rod is governed by the partialdifferential equation. The vibrations of an elastic string are governed by the partial differential equation. However there are complex physical processes which can be described by single system of partialdifferential equations in two or three space dimensions for which it is not always possible to find theclosed form solutions. Often, analytical methods do not exist to solve such problems.In such situations,an alternative is to solve them numerically.Here we study finite difference methods of solving partialdifferential equations.Classification of partial differential equationsConsider general linear partial differential equation of the second order in two independent variables is ofthe form. Such a partial differential equation is said to be2 (a) Elliptic if B - 4AC < 0For example, Laplace’s equation: U+ U= 0 where A = C = 1, B = 0xx yy 2 (b) Parabolic if B - 4AC = 0For example, the heat or di?usion Equation U = ß U where A = 1, B = C = 0txx2 (c) Hyperbolic if B - 4AC > 01 1 For example the 1-D wave equation U =U where A = 1, C = -, B = 0xx tt 2 2 ? ? The classification of the second order linear PDE is very important.It governs the number and nature ofconditions that must be prescribed in the problem in order that a unique solution may exist. Parabolic andhyperbolic problems are classified as initial boundary value problems where as the elliptic problems areclassified as boundary value problems (bvp). Elliptic bvp usually holds inside a closed region (or in anopen region which can be mapped conform ally onto a closed region).Mesh or Grid points.Consider a rectangular region R in x, y plane. Divide this region into a rectangular network of sides 1x =h, Ay = k. The points of intersection of dividing lines are called mesh or nodal or grid points. 5 - 9In difference methods, the partial derivatives in the differential equation are replaced by suitabledifference quotients, converting the differential equation to a difference equation at each nodal point. This procedure is called discretization of the differential equation.The difference equation is applied ateach nodal point.The given data is used to modify the difference equation at the nodes near or on theboundary.The solution of this algebraic system of equations (consisting of the difference equation andthe initial boundary conditions) gives the numerical solution of the given initial or boundary valueproblem.Elliptic EquationsLaplace’s equation?u ? 0 2 2 ? u ? u u(x, y) ? ? ? 0.. 2 2 ?x ?y Poisson’s equation?u ? f 2 2 ? u ? u u(x, y) ? ? ? f (x, y).2 2 ?x ?y Solution of Laplace’s EquationAn alternative to direct solution of the finite difference equations is an iterative numerical solution. Theseiterative methods are often referred to as relaxation methods as an initial guess at the solution is allowedto slowly relax towards the true solution, reducing the errors as it does so. There are a variety ofapproaches with differing complexity and speed. We shall introduce these methods before looking at thebasic mathematics behind them. ? Jacobi’s method? Gauss-SeidelJacobi’s methodThe Jacobi Iteration is the simplest approach. For clarity we consider the special case when dx = dy. Tofind the solution for a two-dimensional Laplace equation simply: 1. Initialize U to some initial guess. ij 2. Apply the boundary conditions. 3. For each internal mesh point set (n) th Letu i , j denote n Iterative value then Jacobi’s Formula is as follow1 (n+1) (n) (n) (n) (n) u =[u + u + u + u ]i , j i-1 , j i+1 , j i , j+1 i , j-1 45 - 10Gauss-SeidelThe Gauss-Seidel Iteration is very similar to the Jacobi Iteration, the only difference being that the newestimate U* is returned to the solution U as soon as it is completed, allowing it to be used immediatelyij ij rather than deferring its use to the next iteration. The advantages of this are: ? Less memory required (there is no need to store U*). ? Faster convergence (although still relatively slow). On the other hand, the method is less amenable to vectorisation as, for a given iteration; the new estimateof one mesh point is dependent on the new estimates for those already scanned. (n) th Letu denote n Iterative value then Gauss-Seidel’s Formula is as followi , j 1 (n+1) (n+1) (n) (n+1) (n) u i , j =[u i-1 , j + u i+1 , j + u i , j+1 + u i , j-1]4 5.2.1Solved ExamplesExample 1: Given the values of u(x, y) on the boundary if the square in the figure. Evaluate the function2 u(x, y) satisfying the laplace equation ? u = 0 by Jacobi method.Solution:First we find out the initial values of u , u , u , u , let u = 0. Then1 2 3 4 41 u = = 1000 ( Daig. formula)1(1000 +0 + 1000 + 2000) 4 1 u2 = = 625 ( Std. formula)(1000 + 500 + 1000 + 0) 4 1 u = = 875 ( Std. formula)3(2000 + 0 + 1000 + 500)4 1 u = = 375 ( Std. formula)4(875 + 0 + 625 + 0)4 1000 1000 1000 1000 2000 u u 5001 2 2000 uu 03 41000 500 0 0(n+1)1 (n) (n) u =1[2000 + u + 1000 + u ]2 3 4 (n+1) 1 (n) (n) u =2[u + 500 + 1000 + u ]1 4 4 (n+1) 1 (n) (n) u =3[2000 + u +u + 500]4 1 4 (n+1) 1 u = (n) (n) 4[ u + 0 + u + 0]3 2 4 First iteration: (put n = 0 in the above results)(1) 1 u = = 11251(2000 + 625 + 1000 + 875) 4 (1) 1 u = = 7192(1000 + 500 + 1000 + 375) 4 (1) 1 u = = 9693(2000 + 375 + 1000 + 500)4 (1) 1 u = = 3754(875 + 0 + 625+ 0)4"