Numerical,Methods,计算方法Lecture13_IVP5_Handout

Lecture6-QuadratureI.HawkeSchool of Mathematics,University of Southampton,UKMA TH3018/6111,Semester1I.Hawke(University of Southampton)Lecture6Semester11/15 Numerical QuadratureThe aim of numerical quadrature is to computebf(x)d xawhere f is a real function of a single variable x.More complex integrals are best evaluated by analytical reduction to a (set of)integral(s)in the form above.Here we will always assume that f(x)isfinite within[a,b],and frequently we will implicitly assume that it is differentiable.The simplest quadratureStandard situation:the value of thefunction is known(or can beevaluated)at a set of points ornodes.The simplest evaluation of thequadrature use the Riemannintegral:a rectangular area iscomputed for each node,and thetotal integral is the sum of theseareas.I.Hawke(University of Southampton)Lecture6Semester14/15 The simplest quadratureStandard situation:the value of thefunction is known(or can beevaluated)at a set of points ornodes.The simplest evaluation of thequadrature use the Riemannintegral:a rectangular area iscomputed for each node,and thetotal integral is the sum of theseareas.improvements change how the function is represented using afinite amount of information.Simplest assumption:know the val-ues f(x j)at nodes x j.This does notgive us complete information aboutthe function.I.Hawke(University of Southampton)Lecture6Semester15/15 Function representationThe simple rule from the Riemann integral is not very accurate.Most improvements change how the function is represented using afinite amount of information.Simplest assumption:know the val-ues f(x j)at nodes x j.This does notgive us complete information aboutthe function.We need information about howthe function behaves away from thenodes.We assume the function be-haves in a suitably simple way,sothat the resulting interpolating func-tion is easy to use.improvements change how the function is represented using afinite amount of information.Simplest assumption:know the val-ues f(x j)at nodes x j.This does notgive us complete information aboutthe function.If the interpolating function is inte-grable then we are done;computethe integrals between each nodeand sum.I.Hawke(University of Southampton)Lecture6Semester15/15 Function representationThe simple rule from the Riemann integral is not very accurate.Most improvements change how the function is represented using afinite amount of information.Simplest assumption:know the val-ues f(x j)at nodes x j.This does notgive us complete information aboutthe function.The choice of interpolating func-tion is definitely not unique.Withquadrature this is not so much of anissue.improvements change how the function is represented using a finite amount of information.A different assumption is that we can compute the function at any point in the interval,but that we wish to do this as little as possible.We then pick the nodes and the in-terpolating function to get the high-est degree of accuracy possible.If the nodes are chosen based on the values of the function f itself then this is an adaptive algorithm.I.Hawke (University of Southampton)Lecture 6Semester 15/15Lecture6Semester15/15Polynomial interpolationThe general class of methods that have the name of Newton-Cotes formulae are simple and robust.They make the following assumptions: The function f(x)is known(or can be evaluated)at afinite set ofnodes{x j},with j=0,...,N,in the interval[a,b].An interpolating function g(x)is found that passes through thesepoints,i.e.g(x j)=f(x j)=f j,j=0,...,N.The interpolating function is a polynomial or piecewise polynomial.The integral is computed from the exact integral of g(x)over the interval.Different methods based on the order of the polynomials used and restrictions on the nodes.Example ConvergenceThe example converges as expected with resolution.I.Hawke(University of Southampton)Lecture6Semester113/15 SummaryWhen working with functions you have to decide how to represent the function using afinite amount of data.If you decide to do quadrature(integration)using thatThe value of the function is known at afinite number of nodes;The interpolating function is a polynomialthen you are using a Newton-Cotes formula.The trapezoidal rule interpolates the function using a linearpolynomial(straight line).The composite trapezoidal rule breaks the interval into many small pieces and applies the trapezoidal rule to each.If each subinterval has the same width h the trapezoidal rule isvery simple and converges∝h2.I.Hawke(University of Southampton)Lecture6Semester115/15。

第二章算法与可计算性理论2.1 算法算法是信息科学的基础理论。

70年代Knuth出版了《The Art of Computer Programming》, 以各种算法研究为主线，确立了算法的重要性。

1974年获得图灵奖。

如同建筑设计、产品设计，总需要一些工程图一样，在生成实际中常需要进行算法设计。

解决计算问题过程：可计算否→算法设计→编程实现→软件系统。

那么什么是算法？很难给出一个完整的、明确的关于算法的定义。

算法（algorithm）这个词最早出现于12世纪，指用数字实现计算的过程，现在是指用计算机解决问题的程序或步骤。

计算方法与算法是密不可分的。

算法一般没有明确的定义，下面几个问题是算法要重点考虑的:●计算速度。

●存储量。

●数值稳定性。

先不给具体定义，只举例说明。

平面区域的三角剖分问题，看其来简单，但实际算法比较复杂。

Delaunay剖分：给定平面上一个点集，存在一种三角剖分，使三角形的最小角达到最大。

算法的定义：是满足下列条件的一系列计算步骤：1.有限性：有限步内必须停止。

2.确定性：每一步都是严格定义和确定的动作。

3.可行性：每一个动作都能够被精确地机械执行。

4.输入：有一个满足给定约束条件的输入。

5.输出：满足给定约束条件的结果。

算法的分析：1.正确性。

一个算法是正确的，如果它对于每一个输入都最终停止,而且产生正确的输出。

调试程序 程序正确性证明。

程序调试只能证明程序有错误，不能证明程序无错误。

2.复杂性。

预测算法对不同输入所需要的资源量。

复杂性测度：时间，空间，I/O。

用途：为求解一个问题选择最佳算法。

时间：输入产生结果需要的原子操作或步骤数。

空间：输入产生结果所需要的存储空间大小。

复杂性可以分：最坏，最好和平均。

3.稳定性。

4.可读性。

如何定性地分析算法的复杂性？这里需要用到下面几个概念：1.同阶函数集。

θ(f (n ))={g (n )|∃c 1,c 2>0,n 0 ,当n ≥ n 0 , c 1f (n ) ≤ g (n ) ≤ c 2f (n )}称为与f (n )同阶的函数集合，注意当n 充分大以后，f (n )是非负的。

NUMERICAL SOLUTION OF PDES 7314.A finite element method for the transport equationWe consider the approximation of the initial value problem:β·∇u +σu =f,in Ω,u =g on Γin (Ω),where βis a unit vector,σ≥c >0andΓin (Ω)={x ∈∂Ω:β·n <0},where n is the unit outward normal to Ω.We shall assume that Ωis a polygon and T h is a triangulation of Ω.Then it is possible to order the triangles {T 1,T 2,···}such that for each k ,the domain of dependence of T k consists of some subset of Γin (Ω)and {T 1,···,T k −1}.With such an ordering,one can develop an finite element method in an explicit fashion,triangle by triangle.Letting P n (T )denote the polynomials of degree ≤n on T ,the Discontinuous Galerkin method for this problem has the form:For T ∈T h ,find u h ∈P n (T )such that u −h =g h on Γin (Ω)and satisfies (β·∇u h +σu h ,v )T − Γin (T )(u +h −u −h )v β·n ds =(f,v )T ,v ∈P n (T ),where for x ∈Γin (T ),u ±h (x )=lim ǫ→0u h (x ±ǫβ).If we solve these equations using theordering discussed above,then u −h is known at the time it is needed to compute the solutionon triangle T .On each triangle,we need to solve a simple square linear system of equations.Note that the solution produced is a piecewise polynomial,but one that is discontinuous across triangle edges.The key to the analysis of this method is the following identity.Lemma 9.Assume that βis a constant vector.Then(β·∇u,u )T − Γin (T )(u +−u −)u +β·n ds=12 Γout (T )(u −)2|β·n |ds +12 Γin (T )(u +−u −)2|β·n |ds −12 Γin (T )(u −)2|β·n |ds.One important implication of this identity is that it easily follows that the linear system on each triangle has a unique solution for σ>0.To see this,we need only show that if f =0and u −h =0,then u h =0.Choosing v =u +h and using the above identity,we get 12 Γout (T )(u −)2|β·n |ds +12 Γin (T )(u +−u −)2|β·n |ds +σ u h 2T =0,and so u h =0.In analyzing this problem,it is helpful to think of u h as evolving in layers S i ,defined byS 0=∅,S i ={T ∈T h :Γin (T )⊂Γin (Ω−∪j<i S j },j =1,2,···.Within a layer,u h can be developed in parallel.We can also define a sequence of fronts F i ,to which u h has advanced after it has been computed in Ωi =∪j ≤i S j .74NUMERICAL SOLUTION OF PDESIn the case when f =0,we also have a very simple stability analysis that can be expressed in the above terms.Theorem 16.12|u −h |2F i +σ u h 2Ωi ≤12|u −h |2Γin (Ω).Proof.Applying our identity with f =0,we obtain 12 Γout (T )(u −)2|β·n |ds +12 Γin (T )(u +−u −)2|β·n |ds +σ u h 2T =12 Γin (T )(u −)2|β·n |ds.Summing over all the triangles in the layer S i and omitting the positive jump terms,we get 12|u −h |2F i +σ T ∈S iu h 2T ≤12|u −h |2F i −1.The theorem follows by iterating this inequality.Also using this key identity,we are able to show thatu −u h L 2(Ω)≤Ch n +1/2 u n +1.Note that this is not an optimal order error estimate,since the best approximation by polynomials of degree ≤n ,would be O (h n +1).15.The finite volume method for elliptic problemsWe consider the approximation of the problem−div(a ∇u )=f in Ω,u =0on ∂Ω.Finite volume methods are based on an approximation of the balance equation − ∂b a ∇u ·n ds = bf dx,valid for any subdomain b ⊂Ω,where n denotes the unit outward normal to the boundary of b .Note this equation can be obtained by integrating the partial differential equation over the subdoiman b and applying the divergence theorem.One approach to the finite volume method,and the one we will discuss,uses a finite element partition of Ω,where the approximate solution space consists of piecewise linear functions,a collection of vertex-centered control volumes,and a test space consisting of piecewise constant functions over the control volumes.More precisely,we begin with a family of triangulations {T h }of the domain Ωand letX h ={v ∈H 10(Ω):v |T ∈P 1,∀T ∈T h }.To construct the control volumes,we let z T denote the barycenter of T and connect z T to the midpoints of the edges of T with line segments.This partitions each triangle T into three quadrilaterals.We use the notation K z,T to denote the quadrilateral in T which shares the vertex z of the triangle T .To each vertex z of the triangulation,we associate a controlNUMERICAL SOLUTION OF PDES75 volume b z consisting of the union of the quadrilaterals K z,T,where the union is taken over all triangles T containing the vertex z.Thefinite volume method is then tofind u h∈X h such that− ∂b z a∇u h·n ds= b z f dx,for all z∈Z0h,where Z0h denotes the set of interior vertices of the mesh T h.There is a large literature on the use offinite volume methods to approximate partial differential equations.An error analysis for the method presented here can be found in the paper“Error estimates for afinite volume element method for elliptic pdes in nonconvex polygonal domains,”by P.Chatzipantelidis and zarov,SIAM J.Numer.Anal.,(42), no5.,2005,pp.1932-1958.。

第一章 误差§1.误差的来源 实际问题——➠建立数学模型—➠确定数值计算方法——➠编制程序上机算出结果模型误差 截断误差或方法误差 舍入误差§2. 绝对误差、相对误差与有效数字(1) 绝对误差与绝对误差限定义： 绝对误差 x x x e e −==***)( .近似值------↑ ↑------精确值通常，由于x 不知道，所以无法得*e ，故估计*e 的上界*ε，即***||||ε≤−=x x e 或 **ε±=x x .↑------称为近似值*x 的绝对误差限，简称误差限。

(2) 相对误差与相对误差限110 ,210021±=±=x x定义： 相对误差 .)(****x x x x e x e e rr −=== 由于x 未知，所以***x e e r ≈； Q **2*****1)(x e x e x e x e −=−，当||**x e 较小时，***x e x e −是**x e 的平方级，可以忽略不计，∴ 取***x e e r=. 与绝对误差类似，只能估计相对误差绝对值的某个上界*r ε，即**||rr e ε≤ ↑------近似值*x 的相对误差限，得（差）。

（好），%10101|)(| %21002|)(|2*1*=≤=≤x e x e r r .(3) 有效数字若近似值*x 的误差不超过某位数字的半个单位，而从该位数字到*x 最左边的那个非零数字（即自左向右看，第一个出现的非零数字）共有n 位，那么这n 位数字都称有效数字，并称*x 具有n 位有效数字。

X XX x L L =*自左向右看，第一个非零数----↑ ↑-----误差不超过该位数的半个单位 例：L 14159.3==πx ，若取近似值14.3*≈x ，则01.0210015.0|)(|*×≤=L x e ，故*x 具有三位有效数字。

(4) 有效数字、绝对误差、相对误差之间关系如何呢？一般(*) )1010(10)1(121*−−−×++×+×±=n n m a a a x L 01≠a ，即n a a a ~ ;9~1:21是.9~0 且1)1(*1021101021||+−−−×=××≤−n m n m x x m m a x a 10)1(||101*1×+≤≤×Q111121***10211010||||||+−+−×=××≤−=∴n m n m r a a x x x e 定理1：若用）（*式表示的近似值*x 具有n 位有效数字，则其相对误差满足不等式 11*1021||+−×≤n r a e 其中1a 为*x 的第一个非零数字。

Lecture 2: Numerical Methods for Differentiations and IntegrationsAs we have discussed in Lecture 1 that numerical simulation is a set of carefully planed numerical schemes to solve initial value problems numerically. Let us consider the following three types of differential equations.dy (t )dt=f (t ) (2.1) dy (t )dt=f (y ,t )(2.2) ∂y (x ,t )∂t =f (t ,y ,∂y ∂x ,∂2y∂x2,...,ydx ,...∫)(2.3)The numerical methods for time integration of these equations will be discussed in the next section (Section 4). Before we conduct the time integration, we need to determine thedifferentiations∂y ∂x ,∂2y∂x 2,... and integrations ydx ,...∫ on the left hand side of the equation(2.3) at each grid point.In this section, we are going to discuss the following four types of numerical methods, which are commonly used in spatial differentiations and integrations.1. Finite Differences (based on Taylor’s expansion)2. FFT (Fast Fourier Transform)3. Cubic Spline2.1. Finite DifferencesFor convenience, we shall use the following notation in the rest of this lecture notes.f ijk n =f (x =i Δx ,y =j Δy ,z =k Δz ,t =n Δt )=f (x i ,y j ,z k ,t n )Given a tabulate functioni :12...N x i :x 1x 2...x N f i :f 1f 2...f NUsing finite difference method, we can obtain derivatives of a tabulated function f . Table 2.1 lists examples of the first-order and second-order finite-difference expressions of[d f /dx ]x =x i , [d 2f /dx 2]x =x i , and [d 3f /dx 3]x =x i . Table 2.2 lists examples of the finite∫,difference expressions of y=f dxTable 2.1. The numerical differentiations based on finite difference methodExercise 2.1.(a) Show that the Central Difference shown in Table 2.1 is a second-order scheme(b) Show that the Forward Difference and the backward difference shown in Table 2.1 arefirst-order schemes., and(c) Determine the forth order central difference expressions of [d f/dx]x=x i, based on the Taylor expansions of the function f.[d2f/dx2]x=x iExercise e the first-order, the second-order, and the forth-order finite-differences expressions to determine the d f/dx, and d2f/dx2, an analytical function f with a fixed Δx. Determine the numerical errors in your results. Compare the numerical errors obtained from different finite differences expressions.Figure 2.1. A summary diagram of the central difference, forward difference, and the backward difference schemes.Table 2.2. The spatial integrations based on finite difference methodExercise 2.3.Use the first order, the second order, and the forth order integration expressions listed in Table 2.2 to determine y (x =π/4) with dy (x )/dx =cos(x ) and boundary conditiony (x =0)=0. Determine the numerical errors in your results. Compare the numericalerrors obtained from different integration expressions.2.2. FFT (Fast Fourier Transform)A function can be expanded by a complete set of sine and cosine functions. In the Fast Fourier Transform, the sine and cosine tables are calculated in advance to save the CPU time of the simulation.For a periodic function f , one can use FFT to determine its differentiations and integrations, i.e.,d fdx=FFT −1{ik [FFT (f )]} fdx ∫=FFT −1{1ik [FFT (f )]} for k >0.Exercise 2.4.Use an FFT subroutine to determine the first derivatives of a periodic analytical functionf . Determine the numerical errors in your results.Exercise 2.5.Use an FFT subroutine to determine the first derivatives of a non-periodic analytical function f . Determine the numerical errors in your results.2.3. Cubic SplineA tabulate function can be fitted by a set of piece-wise continuous functions, in which the first and the second derivatives of the fitting functions are continuous at each grid point. One need to solve a tri-diagonal matrix to determine the piece-wise continuous cubic spline functions. The inversion of the tri-diagonal matrix depends only on the position of grid points. Thus, for simulations with fixed grid points, one can evaluate the inversion of the tri-diagonal matrix in advance to save the CPU time of the simulation.For a non-periodic function f, it is good to use the cubic spline method to determine its differentiations and integrations at each grid point.Results of differentiations obtained from the cubic spline show the same order of accuracy as the results obtained from the forth order finite differences scheme.Exercise 2.6.Use a Cubic Spline subroutine to determine the first derivatives of an analytical functionf. Determine the numerical errors in your results.ReferencesHildebrand, F. B., Advanced Calculus for Applications, 2nd edition,Prentice-Hall, Inc., Englewood, Cliffs, New Jersey, 1976.Press, W. H., B. P. Flannery, S. A. Teukolsky, and W. T. Vetterling, Numerical Recipes (in C or in FORTRAN and Pascal), Cambridg e University Press, Cambridge, 1988.System/360 Scientific Subroutine Package Version III, Programmer’s Manual,5th edition, IBM, New York, 1970.。