CAAM452Lecture10

nullNumerical Methods for Partial Differential Equations Numerical Methods for Partial Differential Equations CAAM 452 Spring 2005 Lecture 9 Instructor: Tim WarburtonTodayTodayIntroduction to a very basic finite volume method Mass ConservationMass Conservation For the following we are going to consider a 1-dimensional domain which we parameterize with the variable x. Now imagine that this line is a figurative representation of a pipe which contains fluid. At every point on the line the fluid has a density measured in mass per meter ( with units kgm-1 ) We define a new function which is a non-negative, real valued function defined on the space-time domain. xDerivation of Mass Conservation LawDerivation of Mass Conservation LawNext we consider an arbitrary section of the pipe, say [a,b] We now assume that the fluid is not created or destroyed at any point inside the section and is traveling with velocity u (which is a function of space and time). For the moment we will assume that u is positive (i.e. the fluid is flowing in the direction of positive x) This allows us to state the following: The time rate of change of the total fluid inside the section [a,b] changes only due to the flux of fluid into and out of the pipe at the ends x=a and x=b. A simple formula relating these two quantities is: nullIn detail: The time rate of change of total mass in the section of pipe [a,b]The flux out of the section at the right end of the section of pipe per unit timeThe flux into the section at the left end of the section of pipe per unit timeUse the Fundamental Theorem of CalculusUse the Fundamental Theorem of CalculusLook carefully at the right hand side: Clearly we may rewrite this as: From which we may deduce:Finally…Finally…Assuming that the integrand of: is continuous and noting that this relation holds for all choices of a,b then we may deduce: In short hand:Advection EquationAdvection EquationLet’s choose a simple, constant, fluid velocity Then the pde reduces to the advection equation: This is a pretty easy equation to solve . Consider the change of variables:nullNote that the units of each variable are consistent. Basic calculus: From which we obtain: Solution and InterpretationSolution and InterpretationSo we know: Which we can instantly solve: where: So an interesting property of the advection equation is the way that the profile of the solution does not change shape but it does shift in the positive x direction with constant velocitySpace Time DiagramSpace Time DiagramLet’s track the information: The dashed lines are which are known as characteristics of the equation. If we choose a point on one of these dashed lines and track back down to t=0 and we will find the value of the density which applies at all points on the dashed line xRecallRecallAssuming that the integrand of: is continuous and noting that this relation holds for all choices of a,b then we may deduce: Well – this does not hold if the density is discontinuous and the integral equation is the appropriate representation:Building a Finite Volume SolverBuilding a Finite Volume SolverLet’s consider the advection equation: Next we take a finite portion of the real line from x1 to xN divided into N-1 equal length sections In each section we will approximate the density by a constant valuex1xNdxPiecewise Constant ApproximationPiecewise Constant ApproximationIn the n’th section the density will be approximated by the constant: x1xNnullChoosing a = xi, b=xi+1, Is approximated (to first order in time) by: where the time axis has been divided into sections of length dt and the i’th cell average at time n*dt is represented by Outstanding question: Now we have to figure out how to evaluate the density at the interval end points given the cell averages. Upwind Treatment for Flux TermsUpwind Treatment for Flux Terms Recall that the solution shifts from left to right as time increases. Idea: use the upwind values Basic Upwind Finite Volume MethodBasic Upwind Finite Volume MethodsimplifyNote we must supply a value for the left most average at each time step:Basic Upwind Finite Volume MethodBasic Upwind Finite Volume MethodApproximate time derivative and look for solution Note we must supply a value for the left most average at each time step:RecallApproximate fluxes with upwind fluxConvergenceConvergenceWe have constructed a physically reasonable numerical scheme to approximate the advection equation. However, we need to do some extra analysis to determine how good at approximating the true PDE the discrete scheme is. Let us suppose that the i’th subinterval cell average of the actual solution to the PDE at time T=n*dt is denoted by Error EquationError EquationThe goal is to estimate the difference of the exact solution and the numerically obtained solution at some time T=n*dt. So we are interested in the error: For the given finite volume scheme dt and dx will be related in a fixed manner (i.e. dt = Cdx for some C, independent of dx). Suppose we let and then the scheme is said to be of order s.Norms and DefinitionsNorms and DefinitionsWe define the discrete p-norms: We say that the scheme is convergent at time T in the norm ||.|| if: It is said to be accurate of order s if:Local Truncation ErrorLocal Truncation ErrorSuppose at the beginning of a time step we actually have the exact solution -- one question we can ask is how large is the error committed in the evaluation of the approximate solution at the end of the time step. i.e. choose Then We expand about xi,tn with Taylor series:Estimating Truncation ErrorEstimating Truncation ErrorInserting the formulas for the expanded q’s:Estimating Truncation ErrorEstimating Truncation ErrorRemoving canceling terms:Estimating Truncation ErrorEstimating Truncation ErrorSimplifying:Final FormFinal FormUsing the definition of qUsing that:Interpretation of ConsistencyInterpretation of ConsistencySo the truncation error is O(dx) under the assumption that dt/dx is a constant.. This essentially implies that the numerical solution diverges from the actual solution by an error of O(dx) every time step. If we assume that the solution q is smooth enough then the truncation error converges to zero with decreasing dx. This property is known as consistency. Error EquationError EquationWe define the error variable: We next define the numerical iterator N: Then: So the new error consists of the action of the numerical scheme on the previous error and the error commited in the approximation of the derivatives.Abstract SchemeAbstract SchemeWithout considering the specific construction of the scheme suppose that the numerical N operator satisfies: i.e. N is a contraction operator in some norm then…Estimating Error in Terms of Initial Error and Cumulative Truncation ErrorEstimating Error in Terms of Initial Error and Cumulative Truncation ErrorError at a Time T (independent of dx,dt)Error at a Time T (independent of dx,dt) If the method is consistent (and the actual solution is smooth enough) then: As dx -> 0 the initial error ->0 and consequently the numerical error at time T tends to zero with decreasing dx (and dt). Specific Case: Stability and Consistency for the Upwind Finite Volume SchemeSpecific Case: Stability and Consistency for the Upwind Finite Volume Scheme We already proved that the upwind FV scheme is consistent. We still need to prove stability of. Stability in the Discrete 1-normStability in the Discrete 1-normSo here’s the interesting story. In the case of a zero boundary condition then we automatically observe that the operator is a contraction operator.Boundary ConditionBoundary ConditionSuppose are two numerical solutions with Then: i.e. if we are spot on with the left boundary condition the N iterator is indeed a contraction.Relaxation on Stability ConditionRelaxation on Stability ConditionPrevious contraction condition on the numerical iterator N A less stringent condition is: Where alpha is a constant independent of dt as dt->0Relaxation on Stability ConditionRelaxation on Stability ConditionIn this case the stability analysis yields:InterpretationInterpretation Relaxing the stability yields a possible exponential growth – but this growth is independent of T so if we reduce dt (and dx) then the error will decay to zero for fixed T. Homework 2Homework 2Due on Friday Recall report presentation policy Time will be given in class on Wednesday to go over questions.Homework 2Homework 2Q1) Solve the pde analytically on the domain Q2) Solve the pde analytically. Q3) Solve the pde analytically and explain what happens to the density along the characteristics.PTOHomework 2 contHomework 2 contQ4) Implement the numerical approximation of: By: Homework 2 contHomework 2 contQ4 cont) For N=10,40,160,320,640,1280 run to t=100, with: dx = (8/(N-1)) dt = dx/6 On the same graph, plot t on the horizontal axis and error on the vertical axis. The graph should consist of a sequence of 6 curves – one for each choice of dx. Comment on the curves. NOTE: For the purposes of this test we define error as: Lecture 5Lecture 5We will define stability for a numerical scheme and investigate stability for the upwind scheme. We will compare this scheme with a finite difference scheme. We will consider alternative ways to approximate the flux functions. Homework 2 correctionHomework 2 correctionQ4) Implement the numerical approximation of: By: Homework 2 contHomework 2 contQ4 cont) For N=10,40,160,320,640,1280 run to t=10, with: dx = (8/(N-1)) dt = dx/6 Use On the same graph, plot t on the horizontal axis and error on the vertical axis. The graph should consist of a sequence of 6 curves – one for each choice of dx. Comment on the curves. NOTE: For the purposes of this test we define error as: