Create Presentation
Download Presentation

Download Presentation
## Numerical methods to solve parabolic PDEs

- - - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - - -

**Numerical methods to solve**parabolic PDEs**Mathematical models: 5° Classification**Classification based on the type of the solution (except, usually, the empirical models) ANALYTICAL SOLUTIONS MATHEMATICAL ANALYSIS ALGEBRA NUMERICAL SOLUTIONS NUMERICAL COMPUTATION • Ex.: numerical methods for PDE • FINITE DIFFERENCE METHOD(FDM) • FINITE ELEMENT METHOD (FEM) • METHODS OF LOCATION**Partial differential equation (PDE)**Second order Partial differential equation: where u=f(x,y,t), a function that satisfies this equation and that is at least two times differentiable, is said to be a solution of the equation. To obtain a unique solution, it is necessary supply some additional information, namely initial (IC) and boundary (BC) conditions. BC must to be of the order N-1 if N is the order of the equation.**Partial differential equation (PDE)**Classical classification for 2nd order PDE’s if = b2 - 4ac < 0 : Elliptic = 0 : Parabolic > 0 : Hyperbolic**Partial differential equation (PDE)**Classical classification for 2nd order PDE’s Laplace’s equation for heat conduction a=1, b=0, c=1 = b2 - 4ac = -4 Heat conduction equation is elliptic.**Partial differential equation (PDE)**Classical classification for 2nd order PDE’s Molar Diffusion a=D, b=0, c=0 • = b2 - 4ac = 02 - 4D 0 = 0 Diffusion equation is parabolic.**Partial differential equation (PDE)**Classical classification for 2nd order PDE’s Convection in conservation laws Firstorder /tor/h a=1, b=v, c=0 a=v, b=1, c=0 = b2 - 4ac = v2 = b2 - 4ac = 1 Convection equation is hyperbolic**Second order Partial differential equation**Classification: Note: An equation can belong to a class in a certain field and another in a different field**Parabolic PDE’s**Parabolic (or diffusion) PDE with one spatial independent variable 2° order PDE 2 independent variables Linear with constant coefficients is a diffusion coefficient t Є [0, +∞] x Є [xA, xB] I.C.: u(t=0,x)=u*(x) all x B.C.1: u(t, x=xA)=uA(t) all t (simple case) B.C.2: u(t, x=xB)=uB(t) all t (simple case) Analytical solution:u=u(t,x) where u is a continuous function of t and x**Analytical vs numerical solution**t 0 xA xB x Analytical solution The region described by these two independent, continuous variables (t and x) is a part of a plane, for the problem under consideration: the length variable, x, varies between xA and xB, and the time variable, t, increase without limit from 0. Numerical solution To obtain a numerical solution , one replaces these continuous variables with discrete variables. When these two continuous independent variables are replaced by discrete variables (also called x and t), the new variables are defined at specific points located in the same region as shown in Figure. The relations between these discrete variables are finite differential equations, and it is these finite differential equations which are solved numerically on a digital computer. u=u(t, x) The index i indicates the position along the x-axis, and the index n is used for the t-axis u=u(tn, xi)**General procedure for solving Parabolic PDE equations**1 0 1,2,… i-1, i, i+1,… N-1, N x=xA x=xB The dependent variable, u, is now a function of two discrete independent variables, x and t. It is therefore necessary to use two subscripts to specify the value of u at a given point; thus: u(xi, tn)=uin**General procedure for solving Parabolic PDE equations**1 0 1,2,… i-1, i, i+1,… N-1, N x=xA x=xB u(xi, tn)=uin The value of the dependent variable is unknown at a row of points at each time level, and there are actually an unlimited number of time levels. It is not feasible to solve for all the unknown values of u simultaneously even when a limited number of time level are considered. Consequently the technique employed is to solve for the unknown values of u at one time level, using the know values of u at the previous time level.**General procedure for solving Parabolic PDE equations**1 0 1,2,… i-1, i, i+1,… N-1, N x=xA x=xB The values of u at the time level when n=0, are given by the initial conditions of equation (IC). These value are used to determine the unknown values of u at the next time level for which n=1. The same computational values is then used to find the values of u for n=2 from the now known values of u at n=1, and so on. I.C. I.C.**General procedure for solving Parabolic PDE equations**t (i,n+1) (i-1,n+1) (i+1,n+1) n+1 n (i,n) x i-1 i i+1 This procedure is continued for as many time increments as desired. Therefore, the finite difference equations are formulated so that they contain values of u at two consecutive levels. The index n is used to designate the last time level at which the value of u are known, and the index n+1 is used to designate the next time level at which variables of the dependent variable are unknown. unknown value unless it is given by a boundary condition. known value**Finite Difference Methods**The finite difference method (FDM) was first developed by A. Thom* in the 1920s under the title “The method of square” to solve nonlinear hydrodynamic equations. *A. Thom and C. J. Apelt, Field Computations in Engineering and Physics. London: D. Van Nostrand, 1961. “The finite difference techniques are based upon the approximations that permit replacing differential equations by finite difference equations. These finite difference approximations are algebraic in form, and the solutions are related to grid points”.**Finite Difference Methods**• Thus, a finite difference solution • basically involves three steps: • Dividing the solution into grids of nodes • Approximating the given differential equation by finite difference equivalence that relates the solutions to grid points (In finite difference methods, each derivative of the PDE is replaced by an equivalent finite difference approximation) • Solving the finite difference equations subject to the prescribed boundary conditions and/or initial conditions**Finite Difference Methods**0 L t (i,n+1) (i-1,n+1) (i+1,n+1) n+1 n (i,n) x i-1 i i+1 If xA=0 and xB=L, the region between 0 and L along the x-axis is divided into N equal increment of size Dx, whit grid points being placed on each boundary. The time-axis is divided into increments of size Dt. tn=n*Dt xi=xA+i*Dx=i*Dx Some useful relations between values of the independent variable at adjacent points: xi+1=xi+i*Dx xi-1=xi-i*Dx For many numerical solutions, it will be desirable to increase the size of the time step, Dt, as the solutions progresses.**Finite Difference Methods**t (i-1,n+1) (i,n+1) (i+1,n+1) n+1 n (i,n) x i-1 i i+1 Both sides must to be evaluated at the same conditions (xi and tn) We will indicate: Evaluate at time tn Evaluate at length xi All i and n The FDM write both sides and check the equality for each grid points**Finite Difference Methods**• The basis of a finite difference method is the Taylor series expansion of a function. • •Consider a continuous function f(x). Its value at neighboring points can be expressed in terms of a Taylor series as: • f(xi+ ∆x) =f(xi) +(df/dx)xi (∆x) + (d2f/dx2)xi (∆x2)/ 2! +.. +(dnf/dxn)xi (∆xn)/n! +…. • •The above series converges if ∆x is small and f(x) is differentiable • •For a converging series, successive terms are progressively smaller**Finite Difference Methods**The terms in the Taylor series expansion can be rearranged to give (df/dx)xi =[f(xi+ ∆x) -f(xi)] /∆x-(d2f/dx2)xi (∆x)/2! -…-(dnf/dxn)xi (∆xn-1)/n! -... Or (2)(df/dx)xi ≈[f(xi+ ∆x) -f(xi)] /∆x+O(∆x) •Here O(∆x) implies that the leading term in the neglected terms are of the order of ∆x, i.e., the error in the approximation reduces by a factor of 2 if ∆x is halved. •Equation (2) is therefore a first order-accurate approximation for the first derivative.**Finite Difference Methods**Other Approximations for a First Derivative Other approximations are also possible. Writing the Taylor series expansion for f(xi- ∆x), we have (3)f(xi-∆x) =f(xi) – (df/dx)xi (∆x) +(d2f/dx2)xi (∆x2)/ 2! -.. +(dnf/dxn)xi (∆xn)/n! + •Equation (3) can be rearranged to give another first order approximation : (4)(df/dx)xi ≈[f(xi) -f(xi-∆x)] /∆x+O(∆x) •Subtracting (3) from (1) gives a second order approximation for df/dx: (5)(df/dx)xi ≈[f(xi+ ∆x) -f(xi-∆x)] /(2∆x) +O(∆x2)**Finite Difference Methods**In summary: • Forward-difference formula (first order-accurate approximation): (df/dx)xi ≈[f(xi+ ∆x) -f(xi)] /∆x+O(∆x) •Backward-difference formula (first order-accurateapproximation): (df/dx)xi ≈[f(xi) -f(xi-∆x)] /∆x+O(∆x) •Central-difference formula (second order-accurate approximation): (df/dx)xi ≈[f(xi+ ∆x) -f(xi-∆x)] /(2∆x) +O(∆x2)**Finite Difference Methods**Geometrical interpretation of finite differences Given a function f(x) shown in Fig. 2, we can approximate its derivative, slope or tangent at P by the slope of the arcs PB, PA, or AB, for obtaining the forward difference, backward-difference, and central-difference formulas respectively. In general, a second order correct analog is always better approximation than is a first order correct analog**Finite Difference Methods**Higher derivates Upon adding (1) and (3) (1) f(xi+ ∆x) =f(xi) +(df/dx)xi (∆x) + (d2f/dx2)xi (∆x2)/ 2! +.. +(dnf/dxn)xi (∆xn)/n! +…. (3) f(xi-∆x) =f(xi) – (df/dx)xi (∆x) +(d2f/dx2)xi (∆x2)/ 2! -.. +(dnf/dxn)xi (∆xn)/n! + we obtain: f(xi+ ∆x) - f(xi-∆x) =2f(xi) + (d2f/dx2)xi (∆x2)/ 2! +O(∆x4) and we have: (d2f/dx2)xi ≈[f(xi+ ∆x) - 2f(xi) + f(xi-∆x)] /(∆x2) +O(∆x2) (second order-accurate approximation) Higher order finite difference approximations can be obtained by taking more terms in Taylor series expansion.**Type of Errors**1)Discretization or truncation error:εi=yi-y(xi) The discretization error encountered in integrating a differential equation across one step is sometimes called local truncation error. This error is determined solely by the particular numerical solution procedure selected, that is, by the nature of the approximations present in the method; this type of error is independent of computing equipment characteristics. Example: when one replace the first derivate by the finite Forward-difference formula,the truncation error is of the order Dx, while it is of the order Dx2 when a Central-difference formula is used.**Finite Difference Methods**Example: • Forward-difference formula (the truncation error is of the order of Dx): (df/dx)xi ≈[f(xi+ ∆x) -f(xi)] /∆x+O(∆x) •Backward-difference formula (the truncation error is of the order of Dx): (df/dx)xi ≈[f(xi) -f(xi-∆x)] /∆x+O(∆x) •Central-difference formula (the truncation error is of the order of Dx2): (df/dx)xi ≈[f(xi+ ∆x) -f(xi-∆x)] /(2∆x) +O(∆x2)**Type of Errors**2)The round-off error is determined, for any particular method, by the computing characteristics of the machine which does the calculations due to its finite memory which lead to a an approximation of any irrational number by a “rounded” value.**Type of Errors**Reducing the mesh size could increase accuracy, but the mesh size could not be infinitesimal. Decreasing the truncation error by using a finer mesh may result in increasing the round-off error due to the increased number of arithmetic operations. A point is reached where minimum total error occurs for any particular algorithm using any given word length.**Finite Difference Methods**Analog finite difference equation for each total derivates • Forward-difference formula (first order-accurate approximation): (df/dx)xi ≈[f(xi+ ∆x) -f(xi)] /∆x+O(∆x) •Backward-difference formula (first order-accurateapproximation): (df/dx)xi ≈[f(xi) -f(xi-∆x)] /∆x+O(∆x) •Central-difference formula (second order-accurate approximation): (df/dx)xi ≈[f(xi+ ∆x) -f(xi-∆x)] /(2∆x) +O(∆x2) Second derivate (second order-accurate approximation) (d2f/dx2)xi ≈[f(xi+ ∆x) - 2f(xi) + f(xi-∆x)] /(∆x2) +O(∆x2)**Finite Difference analog to the partial derivates**Partial derivatives of u at the (i,n)th node Forward Backward Central Forward**Finite Difference scheme**Differential equations Estimating derivatives numerically Finite difference equations (algebraic equations)**Finite Differencing of Parabolic PDE’s**Consider a simple example of a parabolic (or diffusion) partial differential equation with one spatial independent variable with a(diffusion coefficient) constant. xЄ[xA,xB]: I.C.: u(x,0)=uo(x) t>0 : B.C.1: u(0,t)=uA(t) B.C.2: u(xN,t)=uB(t)**Finite Differencing of Parabolic PDE’s**Dt 1 0 1,2,… i-1, i, i+1,… N-1, N x=xA x=xB Dx To apply the difference method to find the solution of a function u(x,t), we divide the solution region in x-t plane into equal rectangles or meshes of sides Δx and Δt. x = i•Δx, i=1,2,…,N t = n•Δt . n=1,2,…**Finite Differencing of Parabolic PDE’s**1 0 1,2,… i-1, i, i+1,… N-1, N x=xA x=xB In the development of the analog finite difference equation we write All i and n Both sides must to be evaluated at the same conditions xi and tn. Then we substitute the analog finite difference equation for each derivates.**Finite Difference analog to the partial derivates**Partial derivatives of u at the (i,n)th node Forward Backward Central Forward**Finite Difference scheme**• Depending on analog finite difference equation used one can obtain different finite difference equations analogous to the: • Euler’s Method (explicit or forward) • Laasonen’s Method (implicit o backward) • Crank-Nicholson Method (implicit)**Finite Difference scheme**Euler’s Method (explicit or forward) The forward or explicit difference equation is probably the most well known, although it is the least efficient of all the possible equations which can be used. We use the forward difference formula for the derivative with respective to t and central difference formula with respect to x.**Finite Difference scheme**t n+1 u(x,t) (i,n+1) n (i-1,n) (i,n) (i+1,n) t xi-1 xi xi+1 x Euler’s Method (explicit or forward) Forward Central**Finite Difference scheme**t n+1 (i,n+1) n (i-1,n) (i,n) (i+1,n) xi-1 xi xi+1 Euler’s Method (explicit or forward) Finite difference equation: i=2…N-1 This equation contains only one unknown value, uin+1, and is written explicitly for this unknown. Therefore, this explicit formula can be used to compute u(xi,tn+1) explicitly in terms of u(xi,tn) for i=2..N-1 (internal grid points). The computation of the values of the dependent variable is thus made one point at a time.**Finite Difference scheme**1 0 1,2,… i-1, i, i+1,… N-1, N Euler’s Method (explicit or forward) Finite difference equation: I.C. i=2…N-1 The values of u along the first time row (n=0 or t=0) can be calculated in terms of the boundary and initial conditions, then the values of u along the second row (n=1 or t=Δt) are calculated in terms of the first time row, and so on:**Finite Difference scheme**Euler’s Method (explicit or forward) The explicit formula is simple to implement, and especially easy to use for computing the value of u at each time level, but the its computation is slow. In general, a numerical solution must converge to the solution of the corresponding differential equation when the finite increments , Dx and Dt, are decreased in size. Analysis as shown that a very restrictive relationship between the size of Dx and that of Dt must be satisfied in order for the solution of forward equation to approach that of the differential equation. The restriction required that: One should use very small Dt and large Dx which results to a stable solution but with low accuracy (in order to minimize the truncation error in the x analogs (O (Dx2)), the size of Dx must be small).**Finite Difference scheme**Euler’s Method (explicit or forward) The many transient problems which have boundary conditions independent of time approach a steady state condition. For these problems: The time increment can be continuously increased and made quite large as the solution progresses towards steady state without causing significant truncation error in the time analog. However, for the forward difference equation, the size of Dt must remain on the order of (Dx)2 for the solution to be stable. Thus, a very small value of Dt must be used for stability even when larger value could be used without causing truncation error. A finite difference equation which does not have this restriction is, therefore, a much better one to use as an analog to the differential equation xЄ[xA,xB]:I.C.: u(x,0)=uo(x) t>0 : B.C.1: u(0,t)=uA B.C.2: u(xN,t)=uB**Finite Difference scheme**Laasonen’s Method (implicit o backward) In searching for a new finite difference equation which does not have a restriction on the size of Dt for stability, we might write the finite differences analogs for the derivates at the new or unknown time which is indexed by n+1.**Finite Difference scheme**t u(x,t) (i-1,n+1) (i,n+1) (i+1,n+1) t n+1 n (i,n) x x xi-1 xi xi+1 Laasonen’s Method (implicit o backward) Note:The backward one is the only possible to be implement**Finite Difference scheme**t (i-1,n+1) (i,n+1) (i+1,n+1) n+1 n (i,n) x xi-1 xi xi+1 Laasonen’s Method (implicit o backward)**Finite Difference scheme**Laasonen’s Method (implicit o backward) (1) i=2…N-1 This equation is implicit; that is, it contains three values of the dependent values , u, at the unknown time level, n+1. For N-1 increments in x and the boundary conditions of equation: B.C.1: u(0,t)=uA; B.C.2: u(xN,t)=uB there will be N-2 unknown values of u at each time level, and there are N-2 equations, one corresponding to each points. In summary: N-2 equations (1) for 2<=i<=N-1 B.C.1 for i=1: u1n+1=uA B.C.2 for i=N: uNn+1=uB**Finite Difference scheme**Laasonen’s Method (implicit o backward) The resulting set of equation represents a linear system in which the coefficient matrix is tridiagonal. If N=7, for each time level, one should calculate uin+1 as: B.C1: u1n+1=uA B.C.2: u7n+1=uB The equation can thus be solved by the Thomas method to obtain N-2 values of uin+1. The same method can be used to obtain the values uin+2 from the now known values of uin+1**Finite Difference scheme**Laasonen’s Method (implicit o backward) • There is no restriction on the size of Dt for stability of the backward difference equation. • The size of Dt can therefore be set, independently of the size of Dx, to minimize the truncation error of the Taylor series in time. • All the values of u at the first unknown time level (n=1) will be computed from the initial values (n=0). • In the solution of these equations on a digital computer, it is usual practice to keep the size of Dx constant throughout the calculation.**Finite Difference scheme**Laasonen’s Method (implicit o backward) If the boundary conditions are of the form of equation: B.C.1: u(0,t)=uA; B.C.2: u(xN,t)=uB so that the transient solution approaches a steady state condition, it is usual practice to increase the size of Dtas the solution progresses in order to obtain the solution in a minimum computing time. This is because as the steady-state conditions are approached, the difference in values of u from one time level to the next diminishes if a constant time increment is used. Similarly, as the steady-state conditions are approached, the size of second and the higher time derivates decrease; thus, the same truncation error will result from a larger Dt as the steady-state is approached.**Finite Difference scheme**Laasonen’s Method (implicit o backward) • The backward difference equation is an efficient one, and it is simple to use. • However, it is only first-order correct in time. • It is desirable to find a second-order correct analog to this derivates first order-accurate approximation second order-accurate approximation