Coursehero >> Ohio >> Kent State >> MATH 72202

Course Hero has millions of student submitted documents similar to the one below including study guides, homework solutions, papers, exam answer keys and textbook solutions.

Solutions Numerical of Partial Di erential Equations Jing Li, Kent State University 1 Solving Laplace equation by nite di erence method Let us consider solving the Laplace equation u = uxx + uyy = 0, on a square domain (x, y) [0, 1] [0, 1] with Dirichlet boundary condition. 2 2 Here = x2 + y2 , is the Laplace operator. We discretize the domain by a mesh of n + 2 by n + 2 grid points de ned as the cross points of n + 2 equally spaced horizontal lines and n + 2 equally spaced vertical line. The coordinates of these grids points are (xi , yi ) = (ih, jh), i, j = 0, 1, ..., n, n + 1, where h = 1/ x = 1/ y = 1/(n + 1) in this example. Therefore, any two neighboring grid points in each direction has the distance h, which is called as the mesh size. By replacing the double derivatives in the Laplace equation by the nite di erence, the discrete Laplace equations at all n by n interior grid points give us n2 equations to determine these n2 unknowns ui,j for i, j = 1, 2, ..., n. At each interior point (xi , yi ), for i, j = 1, 2, ..., n, the double derivatives are approximated by their nite di erences uxx = and uyy = u(xi+1 , yi ) 2u(xi , yi ) + u(xi 1 , yi ) + O(h2 ). h2 u(xi , yi+1 ) 2u(xi , yi ) + u(xi , yi 1 ) + O(h2 ). h2 Denote the approximation of the true solution u(xi , yj ) at each point (xi , yj ), by ui,j , i, j = 1, 2, ..., n. De ne the discrete Laplacian for ui,j by h ui,j = ui+1,j 2ui,j + ui 1,j ui,j+1 2ui,j + ui,j 1 + = 0, h2 h2 i, j = 1, 2, ..., n. The original Laplace equation is approximated by n2 discrete Laplace equations on the nodes, (xi , yj ), i, j = 1, 2, ..., n, as ui+1,j + ui 1,j 4ui,j + ui,j+1 + ui,j 1 = 0, i, j = 1, 2, ..., n, h2 where the boundary conditions gives us the values ui,j for these boundary points where i or j equals 0 or n + 1. For example for a mesh of 3 by 3 lines on the square domain, the only unknown is u1,1 for the interior center point, and we can obtain its value by h ui,j = u2,1 + u0,1 4u1,1 + u1,2 + u1,0 = 0, 1 where the right hand side value is given by these boundary conditions. For a mesh of 4 by 4 lines, the unknowns are u1,1 , u2,1 , u1,2 , and u2,2 . The 4 equations are u2,1 + u0,1 4u1,1 + u1,2 + u1,0 = 0, u3,1 + u1,1 4u2,1 + u2,2 + u2,0 = 0, u2,2 + u0,2 4u1,2 + u1,3 + u1,1 = 0, u3,2 + u1,2 4u2,2 + u2,3 + u2,1 = 0, for i = 1, j = 1, for i = 2, j = 1, for i = 1, j = 2, for i = 2, j = 2, . which can be written as the following matrix form u1,1 4 1 1 0 u 1 4 0 1 2,1 1 0 4 1 u1,2 u2,2 0 1 1 4 = u0,1 + u1,0 u3,1 + u2,0 u0,2 + u1,3 u3,2 + u2,3 In general, each single equation in this 5-points formulation is of the form ui+1,j + ui 1,j 4ui,j + ui,j+1 + ui,j 1 = 0, at i, j. 2 Solving Heat equation We consider solving the heat equation ut = c uxx , with initial condition u(0, x) = u0 (x), and boundary conditions u(t, 0) = 0, u(t, 1) = 0, t 0. 0 x 1, 0 x 1, t 0, We know the exact solution of this ODE is, u(x, t) = u0 (x ct), for any t > 0. Here we look at the numerical methods to nd an accurate approximation of the solution. 2.1 Semi-discrete Methods For semi-discrete methods, we only descretize the space derivatives explicitly, while leave the time derivative, which gives us a ODE to solve. First generate space mesh points to discretize the space interval [0, 1], by equally space points 0 = x0 < x1 < x2 < ... < xn < xn+1 = 1. Here xi = i x, i = 0, 1, 2, ..., n + 1, with x = 1/(n + 1). Then at any point xi at a time t, we have the ODE ut (t, xi ) = cuxx (t, xi ), 2 i = 1, 2, ..., n, (1) where we only consider the ODEs at the interior grid points xi , i = 1, 2, ..., n. On the boundary of the interval, u(t, x0 ) = u(t, xn+1 ) = 0 are given by the boundary conditions. The double derivative uxx (t, xi ) was approximated by the central di erence uxx (xi , t) = u(t, xi+1 ) 2u(t, xi ) + u(t, xi 1 ) + O(( x)2 ), ( x)2 i = 1, 2, ..., n. Replacing the space derivative in equation (1) by the central di erence formula, and denoting the approximation of u(t, xi ) by ui (t), i = 1, 2, ..., n, we have the following set of ODES for ui (t), i = 1, 2, ..., n, c ui (t) = (ui+1 (t) 2ui (t) + ui 1 (t)) . (2) ( x)2 In matrix form, we have, for forward di erence u = u1 (t) u2 (t) u3 (t) . . . un (t) 1 c = 0 2 ( x) 2 1 0 2 1 1 2 ... 0 0 . . . ... 0 u1 (t) 0 u2 (t) c u (t) 3 = Au. 0 2 . ( x) . . 0 2 un (t) (3) This system of ODEs for ui (t), i = 1, 2, ..., n, with initial condition ui (0) = u0 (xi ), can be solved by using a ODE solver. This is called the semi-discrete method, where only discretization in the space is implemented. 2.2 Fully discrete Methods The system of ODEs (3) can be further discretized in time, which gives the fully discrete methods of solving the heat equation. Denote time step size by t, and denote each instant time by tk = k t, k = 0, 1, 2, ... . Denote uk the approximation of u(tk , xi ). i Replace the time time derivative in (2) by the forward nite di erence, we have the following equations, for k = 0, 1, 2, ... , c uk+1 uk i i = uk 2uk + uk , i 1 i t ( x)2 i+1 i.e., uk+1 = uk + i i or in matrix form uk+1 = uk + t c Auk , 2 ( x) k = 0, 1, ... . c t uk 2uk + uk , i 1 i i+1 2 ( x) i = 1, 2, ..., n, i = 1, 2, ..., n. where uk represent the vector of the approximation to the solution u(t, x) at time tk . This is in fact the method of solving the ODEs (3) by using the forward Euler method. 3 The eigenvalues of the matrix A are all between 4 and 0. Therefore from the stability if the forward Euler method, it is required that t 2 4 c ( x)2 = ( x)2 . 2c The Backward Euler method can also be used here. Replace the time time derivative in (2) by the backward nite di erence, we have the following equations, for k = 0, 1, 2, ... , c uk+1 uk i i = uk+1 2uk+1 + uk+1 , i 1 i i+1 2 t ( x) i.e., uk+1 = uk + i i or in matrix form uk+1 = uk + t c Auk+1 , ( x)2 k = 0, 1, ... . c t uk+1 2uk+1 + uk+1 , i+1 i i 1 2 ( x) i = 1, 2, ..., n, i = 1, 2, ..., n. As discussed before for the ODE solvers, the backward Euler method is always stable for all time step size t, and the choice of t is only decide by the accuracy of the algorithm. Another alternative of the fully discrete algorithm is the so called CrankNicolson algorithm. Equation (2) is approximated by the Trapezoid rule uk+1 uk c i i = uk+1 2uk+1 + uk+1 + uk 2uk + uk , i i 1 i+1 i i 1 t 2( x)2 i+1 which lead to the iteration I t c c A uk+1 = I + t A uk , 2 2( x) 2( x)2 k = 0, 1, ... . i = 1, 2, ..., n, and is also stable for choice any of t. 3 Solve Advection equation We consider to solve the following advection equation ut + cux = 0, with initial condition u(0, x) = u0 (x), and periodic boundary conditions u(t, a) = u(t, b), t 0. a x b, We know the exact solution of this ODE is, u(t, x) = u0 (x ct), for any t 0. Now we look at the numerical methods. 4 3.1 Semi-discrete Methods First generate space mesh points to discretize the space interval [a, b], by equally space points a = x1 < x2 < ... < xn < xn+1 = b. Here xi = a + (i 1) x, for i = 1, 2, ..., n + 1, with x = (b a)/n. Then at any point xi at a time t, we have the ODE ut (t, xi ) + cux (t, xi ) = 0, i = 1, 2, ..., n, (4) where we only consider the ODEs for xi , i = 1, 2, ..., n, because u(t, xn+1 ) = u(t, x1 ) from the periodic boundary condition. There are several di erent ways to approximate the derivative ux (t, xi ) by nite di erence formulas: forward nite di erence : backward nite di erence : central nite di erence : u(t, xi+1 ) u(t, xi ) + O( x) x u(t, xi ) u(t, xi 1 ) ux (xi , t) = + O( x) x u(t, xi+1 ) u(t, xi 1 ) ux (xi , t) = + O(( x)2 ). 2 x ux (xi , t) = for i = 1, 2, ..., n, where u(t, xn+1 ) = u(t, x1 ) and u(t, x0 ) = u(t, xn ), from the periodic boundary condition. Replacing the time derivatives in equations (4) by those nite di erence formulation, and denoting the approximation of each u(t, xi ) by ui (t), i = 1, 2, ..., n, we have a systems of ODEs for ui (t). for forward : for backward : for central : c (ui+1 (t) ui (t)) , i = 1, 2, ..., n, x c ui (t) = (ui (t) ui 1 (t)) , i = 1, 2, ..., n, x ui+1 (t) ui 1 (t) c , i = 1, 2, ..., n. ui (t) = x 2 2 ui (t) = , In matrix form, we have, for the forward di erence u = u1 (t) u2 (t) u3 (t) . . . un (t) c = x 1 0 0 1 1 0 1 1 0 1 .. 0 0 . 0 0 0 1 u1 (t) u2 (t) u3 (t) . . . un (t) with the real parts of all its eigenvalues are in the interval [0, 2c/ x]. Therefore if c > 0, i.e., when the wave is moving to the positive direction, this ODE is not stable. 5 For for backward di erence u = u1 (t) u2 (t) u3 (t) . . . un (t) c = x 1 0 0 1 1 0 0 1 1 .. 0 0 0 . 1 0 0 1 u1 (t) u2 (t) u3 (t) . . . un (t) , (5) with the real parts of its eigenvalues in the interval [ 2c/ x, 0]. Therefore if c < 0, i.e., when the wave is moving to the negative direction, this ODE is not stable. For the central di erence, u = u1 (t) u2 (t) u3 (t) . . . un (t) c = 2 x 0 1 0 1 1 0 0 1 1 0 ... 0 0 1 1 1 0 0 u1 (t) u2 (t) u3 (t) . . . un (t) , 0 with all pure imaginary eigenvalues. Therefore this ODE is stable, no matter the sign of c. The reason why the sign of the wave speed c determines the stability of the ODE obtained in semi-discrete algorithm can be explained as follows. Let us assume c > 0, i.e, the wave will move to the right. The solution of the advection equations is simply u(t, x) = u0 (x ct). Consider a family of lines x ct = constant = on the two dimensional space (x, t), where each line corresponds to one constant . We call each line as a characteristic line. On each characteristic line, the function value is the same, which means the solution information is in fact propagating along each of characteristic line. Suppose we know the value u(t, x) at the time t, for each i = 1, 2, ..., n, then after an interval of time t, the solution u(t + t, x) is fact determined by the value u(x c t, t) at the previous time t, since the two points (x, t+ t) and (x c t, t) are in the same characteristic line. This means that when we discretize the space derivative in the ODE ut (t, xi ) + cux (t, xi ) = 0, at each xi , it would be better to include some information from the left side of xi . Therefore both backward di erence formulation and central di erence formulation give us stable ODE to solve, since we include in the formula the information from the left point u(t, xi 1 ). While for the forward di erence formulation, which use u(t, xi ) and u(t, xi+1 ), it is unstable. 6 3.2 Fully discrete algorithm In the fully discrete algorithm, time derivatives in (4) are also discretized by nite di erence formulas. Here we assume c > 0 and we consider to solve the semi-discrete ODE ui (t) = c (ui (t) ui 1 (t)) , x i = 1, 2, ..., n, (6) which are derived by using backward nite di erence formulation to discretize the space derivatives and are (5) in matrix form. Denote the time step size by t, and tk = k t, k = 0, 1, ... . Denote by k ui the approximation of u(tk , xi ). Using forward Euler method to solve the system of ODEs (6), we have uk+1 = uk i i c t k u uk , i 1 x i i = 1, 2, ..., n. (7) Since the real parts of the eigenvalues of the matrix in (5) are in the interval [ 2c/ x, 0], we know, from the stability of the forward Euler method for solving ODE, that the step time-step size t need satisfy t x 2 = . 2c/ x c We can also use another approach here to show that the stability requires that t x/c. We rst show that t x/c is su cient to for the stability of forward Euler method (7). Let us denote uk the solution vector (uk , uk , ..., uk ) at the time tk . We de ne the norm of uk by 1 2 n u k n = i=1 uk+1 . i Using the forward Euler formulation (7), we have u k+1 n = i=1 n uk+1 i n = i=1 uk i c t k ui uk i 1 x i=1 c t k c t k ui + u. 1 x x i 1 If t x/c, then we know 1 c t/ x 0, therefore uk+1 1 c t x n i=1 uk + i c t n k u = uk . x i=1 i 1 Therefore if t x/c, then the forward Euler is stable. In the following we will show that if t > x/c, then the error for the forward Euler methods (7) can blow up. Here we give an error propagation example for solving ut + ux = 0, 7 by using the forward Euler method (7). Denote the di erence between uk and i the exact solution value u(tk , xi ) by and assume there are no errors at other place at time tk . From (7), we have uk+1 = uk i i t k ui uk , i 1 x i = 1, 2, ..., n Denote = t/ x, and we have uk+1 = uk (1 )uk . i i 1 i We can see that, if < 1, e.g., = 1/2, then the propagation of the error is in fact decreasing. However, if > 1, for example = 2, the error will become ampli ed. 8

Find millions of documents here - Study Guides, Homework Solutions, Papers, Exam Answer Keys and more. Course Hero has millions of course related materials that will enable you to learn better, faster and get an A in all your courses.
Below is a small sample set of documents: