Randall J.LeVeque——Finite Difference Methods for Ordinary and Partial Differential Euqations Mark H.Holmes——Introduction to Numerical Methods in Differential Equations 魏朝祯——Lecture in UESTC: 微分方程数值解
In this article we ignore the effect of rouding error arsing in float number computing.
Basic definitions
Truncation Error (): error of difference method used to approx derivative;
Local Truncation Error (LTE, ): error of approximation between discrete equation and real equation;
Global Error (): error of using numerical solution to approx real solution;
One-step (global) Error (): error of using numerical solution to approx real solution by a single step;
We use these definitions to discribe the properties of a FDM and its behavior when using in a specific equation in the following. This introduction is divied into two parts according to kinds of BC.
BVP
If the value of solution on the boundary of region considered is given, then the solution of this prob is called a steady solution, and this kind prob is so called boundary value problem (BVP for short). In this case, we only need to select the size of gird in space, which makes it easier compared with problems involved with time.
One of the easiest BVP is that is, a 2-rd ODE with Dirichlet boundary. We take this prob as an example to review the process of using FDM to sovle a BVP:
Select a difference method to approx the derivative;
Replace the derivative in equation by the difference formula in step 1;
Drop the LTE in the discrete equation and replace the value of real solution by grid value;
Combine the difference equations in every grid point and rewrite this system as matrix form;
Sovle the linear system to get a numerical solution on grid. Using the equispaced grid in : and the central difference in step 1: where the truncation error , revealing this difference method has 2-rd arrcury. Using this difference format in step 2 we have the discrete equation: and too, is the LTE. Then we have the discrete equation for computing: now we collect thest equation to have a linear system: where and is the numerical solution and value of on grid respectively.
Obviously, the error arising in step 3 since we drop the LTE in discrete equation to be able to construct a linear system. Denote the values of real solution respect to this problem on the same grid is , then we have and the global error satisfy the following error equation The goal of numerical computation is to have a numerical solution with global error as small as possible in the meaning of some norm. If is nonsingular, then . In this example, has 2-rd accuracy, but it’s possible that has some eigenvalue near , which would expand the LTE, and this eigenvalue may be even closer to 0 as we increase leading to low accuracy of final solution. In fact, it could be better to rewrite the error equation as to emphasize are dependent on grid size . To ensure is small, a natural idea is assuming is small and is not singular as . In this way of thinking, we are leaded to the concepts of consistency, stability and convergence.
Frankly speaking, we call a FDM is consistent if , and stable if absolute value of minimal eigenvalue of is bounded below when , and finally convergent if global error . It turns out that the consistency and stability of FDM is sufficient to guarantee its convergence. This statement is also called fundamental theorem of finite difference method (FT of FDM).
Convergence Analysis
Due to the importance of these concepts, the precise definitions are as followed. A FDM for a BVP is consistent if as ,which means the discretization is sensible to the problem; stable if exists for all sufficiently small (for ) and if there is constant , independent of , such that for all , which assume is not too bad; convergent if global error as , which means numerical solution is a good approximation of real one if is sufficiently small. By this concepts, the fundamental theorem of finite difference method is In fact, we can go further to say in the case of BVP by previous analysis.
Remark
Usually, the consistency of a method is easy to check, however, the stability is much more harder to analysis. The main difficulty of stability analysis is number of eigenvalue increases as increases, and not easy to give exact form of these eigenvalues in general. For some complex problems, it may not even be clear how to define stability in an appropriate way to make statements like FT of FDM hold, in this case, the challenge is to find a definition of “stability” that allows one to prove convergence by things like FT of FDM while at the same time being sufficiently manageable that we can verify it holds for specific FDMs.
IVP
There is a basic difference between BVP and IVP:
BVP—parallel compute in space, give answer in a single step;
IVP—serial compute in time while parallel compute in space, solution of some time depends on before solutions;
This difference makes the error analysis of BVP and IVP different. Generally, error analysis of IVP is more complex than BVP.
In the field of IVP, “stability” means the amplification factors of each one-step, which contribute to the global error, can be bounded uniformly.
ODEs
Problems with first order ODE system and initial value is one of the most common IVP. The IVP in this section takes the form where could be a vector value function. Like before, we need to define convergence in this case first. To discuss the convergence of a FDM for this kind of IVP, we focus on a fixed (but arbitrary) time and consider the error in our approximation to computed with the method using time step . Let us denote the computed value at as , where .
When talking about convergence in IVP, it’s important to make a distinction between convergence of a method and convergence of a method in a specific problem. To speak of a method being convergent in general, we mean that it converges on all problems in a reasonably large class with all reasonable staring values, while in the later case we only need to consider this method applied to the single problem.
Here we consider the convergence of general r-step linear method, which has the following form in general in the brackets indicate the time step is , hence we need staring values to employ a -step LMM (linear multistep method). Now we give the definition of convergence.
An -step LMM is said to be convergent if applying this method to any ODE like before with Lipschitz continuous in (to ensure solution exists), and with any set of starting values satisfying we obtain convergence in the sense of for every fixed time at which the ODE has a unique solution.
Note this definition is quite strict and hard to verify directly in practical, which assume the LMM satisfying the condition for a big class of problems. Fortunately, with the help of ODE theory, to check whether a LMM meets the condition or not, we can apply it to a special class of problems called test problems to verify.
Constant Coefficients Linear OEDs
In the beginning of this section, the ODE system is which may not be a linear system. By using Taylor’s theorem, we consider its first order approximation: It’s sufficient to discuss the constant coefficients linear ODEs to describe the properties of original problem in if is sufficiently small, where could be a matrix if is vector-valued. We know that this ODEs has a special solution: In linear algebra, the generalized eigenvector space decomposition theorem tells us this solution can be wrote as where is eigenvalue of , is the algebra multiple of , and is the corresponding component of . Note that can be expanded by , we are leaded to consider the simple scaler linear equation: with the initial data . This is the test problem we need.
Zero-Stability
In short, a method is zero-stability means it is “stable in the limit as ”, so that it would not amplify the errors arising in previous steps. We take the analysis of one-step methods as an example to put this clear, then give the criterion for the general r-step LMM.
One-Step Method
Use the forward Euler method to approx derivative to solve the test problem we have then the error arising in this discretization satisfies where is the LTE, and we denote as , that is, the one-step error of -th step. Then the global error is If , the first term in RHS increasingly tends to as , while the second term can be bounded by increasingly as . Combine these terms we have If , we have similar inequality. In both cases we have by using and . Note , the global error is bounded by multiple of LTE, thus this method is convergent.
In the global error -th step error is amplified by the factor , the stability of method shows up in the inequality uniformly, that is, the amplified factor is bounded independent of as , which makes us to conclude that each contribution to the final error can be bounded in terms of its original size as a one-step error, hence the global error has the same order of magnitude as the LTE.
r-step LMM
Apply the general r-step LMM to the test problem we have and the error equation is Like our analysis in the case of one-step method, we expect the global error given by a zero-stable method can be bounded in terms of the sum of all the one-step errors and hence has the same asymptotic behavior as the LTE as . Since can be arbitrarily small, we consider the limit situation , and the error equation deduced to If the method is convergent, the solution () of this linear system must be at least bounded, which assumes the solution of homogenous system to be bounded, that is, sequence is a bounded sequence. To solve the sequence from above recurrence relation, there is a general approach by using generating function. Precisely, the generating function is a formal power series: We need to use to construct the a new series whose coefficients of each term happen to coincided with our target by consider where if . Use the recurrence relation we have as a consequence, has a rational representation: where are known initial values. Expand this rational function into power series gives every coefficient by these initial values and . However, we don’t need to really calculate each but only want to check the boundedness of this sequence, so we prefer another method based on matrix.
We construct the following matrix: and denote , then the recurrence relation give use this formal repeatedly we have Let be the Jordan form of , that is, and where is Jordan block , then we have . As a consequence, the first entry of can be written as where are all distinct eigenvalues of and is the algebra multiple of . The coefficients are determined by the initial values given by . To ensure , we assume if and if for each . Since is the companion matrix of polynomial by letting and we have is the characteristic polynomial of , this condition is also called root condition:
for each root of ;
if is a multiple root of , then .
Now we give the definition of zero-stability: A -step LMM is zero-stable if roots of , which is determined by the LMM, satisfy above root condition.
In fact, we can analysis directly instead of and have where and . From this equation, we see that zero-stability guarantee the one-step errors arising in previous steps will not be amplified uncontrollably as , and this ensures the final global error will be bounded by , where is some constant, in this limit situation ().
However, zero-stability only tells us the method would work in the case of , that is, , which is impractical. A zero-stable method may be unstable if .
Absolute Stability
How small should be to guarantee a small global error? Absolute stability of a method helps us to choose an appropriate .
We consider the general -step LMM again and have but now and . Denote , where can be a complex number. We can still construct the matrix but denote it as since now it depends on instead of being a constant matrix like in previous section. And our polynomial turn to which is still the characteristic polynomial of , where then the root condition becomes
for each root of ;
if is a multiple root of , then . The region in complex plane consisted by points satisfying root condition is called the (absolute) stable region of this method. Note a -step LMM is zero-stable if and only if the origin point is located in its stable region.
The error equation in this case is The (absolute) stable region tells us which could be chosen to guarantee the one-step errors arising during the time passing will not be amplified uncontrollably and contribute to the final global error.
There is a useful way to plot the boundary of stable region of a LMM. Let , then gives which is the parametric curve corresponding the boundary of stable region in the plane.
A-stability
There are some problems whose solutions are smooth and slowly varying in most part of the time interval considered, but in a context where the nearby solution curves are much more rapidly varying, which lead to ( first derivative of solution) changing greatly, i.e. varying largely for some fixed . This property force us to use method with as large as possible stable region to solve such problems to ensure stability.
For ODEs, the stiffness of the system is described by the ratio where s are eigenvalues of Jacobin matrix . The bigger is, the more stiff this system is.
To overcome the difficulty caused by stiffness, methods with unbounded stable region are more prefer in numerical. However, it’s proved that all explicit methods can not satisfy this condition, their stable regions are all bounded, we have to use implicit method.
Note if real part of in test problem is negative, then the solution should decrease to fast. It’s natural to assume our method should work well in this case, which means that the stable region of method includes left half plane, and we say such methods are A-stable. To put it clear, an ODE method is A-stable if left half plane is included in its (absolute) stable region. In fact, requiring a method to be A-stable is relatively hard, i.e. it turns out to be an A-stable method can have at most second order accuracy, and the Trapezoidal method is the most accurate method among A-stable methods.
Since the restrictions of A-stable are too hard to satisfy and s of many stiff problems are mainly located near the real axis, the concepts of -stable may be more useful. A method is said to be-stable if its stable region includes .
L-stability
A method with A-stability ensure a small global error with a relatively large , which makes such methods useful to solve stiff problem. But for more difficult problems e.g. problems which are extremely stiff, we need more stable method. In this case, we require the stable region not only includes the left half plane, but also includes the infinity point as an inner point, such methods are defined to be L-stable. For example, the Trapezoidal method is A-table but not L-stable, while Backward Euler method is both A-stable and L-stable.
PDEs
For problems described by PDEs, we need to do discretization in both time and space, which makes the error analysis in this case based on, but very different from the one of ODEs. Here we focus on a class of method for solving these problems called MOL.
Method of Line
One of the most common way to solve PDEs numerical is called method of line (MOL), which deduce the problem of solving PDEs into the problem of solving a system of ODEs by discretize the space firstly. For example, we consider problem given by where is some differential operator about space variable . Generally, MOL first discretize space by using some FDM and have the ODE system: where by using as the space grid points, . The we can use the theory introduced in last section to do error analysis. But note that is a matrix whose size is related with , so the stable region will also change if varies when we use some FDM to discretize time too. This fact means the stability of FDM in the case of PDE not only depends on the property of method, size of and , but also is concerned about the relation between and .
Like in the previous section, we use the LMM to solve this ODEs system and have where is the numerical approximation of and the matrix is dependent on the FDMs we using to discretize space and time. For example, if we use forward Euler method to discretize time, then .
Lax-Richtmyer Stability
Just like the analysis for ODEs, we have where is one-step error arising in -th step. And we get once time again with , we hope the second term can be bounded by some again, and this is what Lax-Richtmyer stability about: A linear method is said to be Lax-Richtmyer stable if , for each time , there is a constant such that for some norm.
If a method satisfies , then it’s obviously Lax-Richtmyer stable, and such methods are called strong stable in fact. It’s easy to check that we can loose this condition to only assume for some constant and have then can serve as of this method.
However, to make use of , we have to figure out all eigenvalues of , which is not easy in general when the problem is not trivial. Fortunately, we have a more practical method to do the stability analysis in some condition.
Von Neumann Analysis
There is a very useful tool for stability analysis called Von Neumann analysis which is based in Fourier Analysis and usually applied to problem with periodic boundary condition or Cauchy problems (problem on the whole space, without boundary).
For now, we consider the Cauchy problem in the whole real axis, and our space grid are , where Our goal is find the region of by which the method we used is Lax-Richtmyer stable.
It’s sufficient to find s satisfying With the help of Parseval’s theorem, this inequality is equivalent with where is the Fourier transform of on the discrete space grid (according to Fourier Transform And Fourier Series), that is which is a scalar (function) instead of a vector in , and the norm deduce to the absolute value (integration on ) of input. What really matters is the scale of this problem has no relation with size of but only depends on the number of PDEs in the original problem. This frees us from the dilemma of finding spectral radius, which is very difficult due to the huge size of if is small.
In practical, Von Neumann analysis focus on a single wave number , this works benefited from the linearity of system. lets (in fact , omit the constant for convenience) and suppose , then take these terms as into the discrete equation given by FDM used to have a representation of , finally request to satisfy with some constant . This will give us the relation between and that they need to satisfy to make the FDM stable. What we do in Von Neumann analysis in fact ensures some inequality like holds by using Cauchy inequality, and this guarantees the stability. There is an example in the MathStackExchange which shows how things work in the case of considering inhomogeneous system.
