PyModPDE：修正方程分析的Python软件

软件X 12（2020）100541原始软件出版物PyModPDE：一个用于修正方程分析的Python软件放大图片作者：James C.托尼·萨阿德·萨瑟兰美国犹他州盐湖城犹他大学化学工程系ar t i cl e i nf o文章历史记录：收到2020年收到修订版2020年5月11日接受2020年保留字：修正方程分析有限差分格式稳定性分析数值方法误差分析a b st ra ct修正方程是分析偏微分方程数值方法的一个有用工具。它提供了一个给定的数值方案的稳定性，扩散和色散特性的洞察。然而，由于涉及大量的代数，PyModPDE是一个python软件，它使用一种新颖的方法来生成修改后的方程式 它以离散偏微分方程为输入，以LATEX格式输出修正后的方程。我们讨论了PyModPDE所基于的新方法，然后使用一维和二维PDE验证软件。PyModPDE是计算科学家的重要工具和工程师，用于教育和研究目的。©2020作者由爱思唯尔公司出版这是CC BY许可下的开放获取文章（http://creativecommons.org/licenses/by/4.0/）中找到。代码元数据当前代码版本v1.0用于此代码版本的代码/存储库的永久链接https://github.com/ElsevierSoftwareX/SOFTX_2020_22Code Ocean compute capsule NA法律代码许可证MIT许可证使用git的代码版本控制系统使用Python的软件代码语言、工具和服务编译要求，操作环境依赖性标准Python 3安装与sympy如果可用，链接到开发人员文档/手册https://github.com/saadgroup/pymodpde问题支持电子邮件tony. utah.edu1. 动机和意义修正后的方程是一个有价值的工具，在分析有限差分近似的时间依赖的偏微分方程（PDE），因为它提供了洞察离散化的影响。虽然适用于线性和非线性偏微分方程[1-这并不奇怪，因为修改后的方程提供了定性和定量的信息，一个给定的有限差分格式的行为和稳定性。修正方程的思想首先由Warming和Hyett [5]在1974年的一篇里程碑式的论文中提出。他们的论点很简单：由于截断误差，*通讯作者。电子邮件地址：mokbel. chemeng.utah.edu（M.卡拉姆），詹姆斯.萨瑟兰@ utah.edu（J.C.萨瑟兰），托尼。萨阿德@ utah.edu（T。Saad）。对于一个给定的微分方程，它并不解那个方程，而是解一个完全不同的偏微分方程。他们的方法包括将有限差分公式中的每一项替换为泰勒级数展开式，然后根据控制方程的规定，将高阶时间导数表示为空间导数。修正后的方程中的项可以用来解释所用有限差分格式的耗散和色散行为（尽管这在非线性偏微分方程的情况下不是很清楚）。尽管它的机械性质，推导出修改后的方程是相当繁琐，并借给自己的用户和程序明错误，由于涉及大量的代数1990年，Chang [6]提出了第一个1替代获得以简单且更系统的方式修改方程，其易于以符号代码实现。他的方法适用于线性时变偏微分方程，1据作者所知。https://doi.org/10.1016/j.softx.2020.1005412352-7110/©2020作者。 由Elsevier B.V.出版。这是一篇开放获取的文章，使用CC BY许可证（http://creativecommons.org/licenses/by/4.0/）。可在ScienceDirect上获得目录列表SoftwareX期刊主页：www.elsevier.com/locate/softx2M. 卡拉姆，J.C.Sutherland和T.Saad/SoftwareX 12（2020）100541=+L==∞∑∑jXy z∑nnu（=X=yz⏐根据冯-诺依曼稳定性分析得出的放大因子与修正方程的系数直接相关目前尚不清楚为什么张修正方程的系数（an）可以通过对方程（1）求微分来获得。（3）关于kn，k=0因未得到人们注意以类似的方式，Ramshaw [7]独立地使用放大因子来估计数值扩散。他1an= inn！αn=0（四）方法包括取放大因子相对于波数k的一阶和二阶导数，在k0处进行评估。Carpentier [8]提出了另一种使用放大因子获得修正方程系数的方法，但可以证明其与Chang [6]提出的方法大致相似在本文中，我们认为Chang我们首先介绍张的工作，在一个简洁的形式，然后将其推广到两个和三个空间维度。然后，我们讨论我们的软件实现的方法，然后给出说明性的例子来验证软件。1.1. 拟议框架对于给定的形式为ut（ux，u xx，. . . ）0的常系数差分格式δtu+γ1δxu+γ2δxxu+···=0，（1）修正方程的形式为这表明，如果α是一个已知量，那么修正方程的系数αn可以通过直接对α求导来计算。我们认为，而不是通过使用泰勒级数和重复微分和次幂来推导修改后的方程，可以改为执行以下操作：1. 替代使用eαtei kx转换为方程中的差分格式。（1）、2. 从等式中求解α（1）、3. 使用等式2找到修改的等式的系数an。（四）、这种方法是直接的，并且可以在Sympy这样的符号包中进行简单的符号操作。1.2. 推广到更高的维度类似的分析可以应用于更高的维度。在这种情况下，我们提出以下形式的修正方程ut =a100ux+a010uy+a001uz+a110uxy+a101uxz+a011uyz+···ut= a1 ux +a2uxx +a3u xxx∞+···=an=1阿努乌nxn（二）=j，∑l，m≥0j+l+mujlmxjylzm.其中ai一般不等于γi。从本质上讲，修正后的方程考虑了由j+l+m≥1（五）数值格式这个方程代表了差分格式实际上代表的物理性质。奇导数代表平流和弥散，而偶导数对应于扩散。基于系数ai，然后可以对有限差分方法的稳定性进行合理的论证，并解释特定数值方案中的数值扩散和色散的存在。修正方程的系数传统上是通过将有限差分公式中的每一项代入泰勒级数而得到的。高阶时间导数，然后通过取代原来的PDE空间导数取代。另一种方法[6，7]是基于这样的事实，注意交叉导数项的存在，例如uxy。这些项在本质上通常是各向异性的，因为它们在不同方向上导致不同的物理行为（在第3.2节中给出了一个例子）。假设一个基本解的形式uk=eαtei（kxx+kyy+kzz）=eαteik·x，（6）我们代入Eq。（5）恢复了放大指数与修正方程∞放大因子与修正方程的系数直接相关。放大因子是从冯-诺依曼稳定性分析中获得的，并指定增长率α=j，l，m≥0+l+m≥1ajlmij+l+mkj kl km.（七）引入到PDE的有限差分表示中的误差，如由等式2给出的。（一）.von-Neumann稳定性方法系数ajlm通过下式获得：p∑通过分析形式u=a=1则该方案是不稳定的，否则是稳定的。在这方面-αjl mk=0（八）keαt eikx，服从等式（一）.如果解是无界的，JLMinj！l！m！∂k∂k∂k⏐然后，放大因子被定义为相对在一个时间步长内，Gu（t+t）即，并指示溶液随时间的放大/衰减。很明显，假设λt>0，放大指数α决定了解随时间增长（α >0，不稳定）还是衰减（α0，稳定）。最后，考虑到PDE是线性的，它足以研究由于叠加的单个频率uk eαt eikx的放大。”[5]“以”“以。（2），可以通过以下公式将放大指数α与修正方程的系数联系起来：∞α=我也是。（三）n=1其中n=j+l+ m。2. 软件描述PyModPDE是用python编写的，因为（1）它在研究人员和学生中很受欢迎，（2）它的易用性和快速的开发时间，以及（3）强大的符号库的可用性，特别是Sympy。 PyModPDE软件提供了一个一阶（时间）时间相关PDE类的单一接口，允许用户构建其右侧（RHS），并以LATEX形式生成修改后的方程表示。一.M. 卡拉姆，J.C.Sutherland和T.Saad/SoftwareX 12（2020）1005413+Δt−= −=⏐x=i、ji− 1，j.（= −DiferentiEquation（dependentVarName，indendentVarsNames，indices =[i，j，k]，timeIndex=n）i，j，n=symbols（DE=DiferntialEquation（depentVarName=2.1. 软件构架PyModPDE 遵 循 的 抽 象 是 基 于 将 PDE 解 释 为 由 左 手 侧（LHS）和右手侧（RHS）组成。LHS总是由因变量φt的一阶时间导数组成。另一方面，RHS由任意数量的具有任意阶的空间项组成，例如aφx或µφxx。为了构造PDE的离散等价物，必须为RHS上的每个项提供离散化方案将LHS的离散化更高阶的离散化是完全可能的，并且将成为代码未来版本的一部分。基于此，PyModPDE定义了一阶线性PDE的类，其中LHS固定为一阶时间导数。这个类允许用户用用户指定的因变量和自变量的名称然后，用户继续创建术语以添加到RHS。当创建项时，用户提供空间导数的阶数、方向（x， y， z）、应该发生离散化的时间级别（对于显式或n1表示隐式）和离散化的类型（稍后讨论）。最后，进行简单的调用以生成修改的方程，并且LATEX格式的方程变得可用。PyModPDE功能的图形描述如图所示。1.一、2.2. 软件功能和示例代码片段使用PyModPDE很简单，需要四个步骤：(1)实例化DifferentialEquation类，（2）构建RHS项，（3）使用在步骤（2）中生成的表达式设置PDE的rhs，并且最后（4）生成修改的方程或放大因子。要讨论API，最好从一个具体的示例开始，之后将详细介绍所有源代码。考虑获得对流方程的修正方程，ut aux在一维空间中具有时间上的前向欧拉，空间上的UPWIND（FTUS方案）离散化。在具有索引i的任意空间网格点处的离散方程-不要与虚单位i混淆-由以下公式给出：将被定义为sympy符号，而PDE的其他属性（如因变量和自变量的名称）则作为字符串传递 在声明所有将要使用的符号（a、i和n）之后，实例化一个DifferentialEquation对象，它表示一阶时间相关微分方程。 此类具有以下构造函数签名哪里dependentVarName是依赖变量名称的必需字符串。independentVarsNames是独立变量名称的必需字符串列表。索引是对应于用于指代不同网格位置的索引的合数符号（例如，[i，j，k]）的列表，例如，ui，j. 这些符号按照指定的顺序与独立变量匹配。timeIndex是时间索引的符号，例如n。请注意，索引和timeIndex参数中使用的符号必须事先声明。作为构造二维量uu（t， x， y）的微分方程的一个例子。我们将使用i和j作为空间索引，分别对应于x和y，n作为时间索引。一旦DifferentialEquation对象被实例化，下一步就是开始在RHS中为这个PDE构建项有两种方法可以实现这一点。从第一个方法开始，可以使用成员函数expr。这个方便的函数自动生成一个有限差分近似的基础上泰勒级数展开的n阶空间导数在一个特定的方向，并返回其离散化的符号表达式。函数expr有n+1−nnnn以下签名伊伊阿Δtii−1.（九）阿克斯在PyModPDE中，可以如下实现哪里order是表示导数的空间阶数的必需整数。directionName是独立变量名的必需字符串，我们正在区分该变量。time是我们评估该表达式的时间级别所需的Sympy符号，例如n+1，n，\ldots。 模板是在网格点i处具有参考索引0的有限差分表示中使用的N个点的所需整数列表。注意，阶数必须小于N，有限差分表示才存在。例如，假设要使用后向差分公式在u（t， x， y）的x方向上添加平流项uun−un这个例子说明了我们前面提到的四个步骤。通过设计，PDE的所有系数以及空间和时间索引要使用expr实现这一点，可以设置order=1（一阶导数）、directionName=fromsrc. pymodpde imp ortDiferentienEquationfrenfre n f symbols#declaretheadvectionvelocitya=symbols（“a“）#declarethespatialanddtemporalindicesi，n=symbols（“in“）#construct ttimeded ededetifferntialeqationDE=DiferntieEqation（ded e ded entVarName =“u“，independentVarsNames=[“x“]，indices =[i]，timeIndex=n）#我是如何处理这些问题的（我将如何处理这些问题 在此区域中的位置）：advectionTerm=-a*DE. expr（order =1，directionName=“x“，time=n，stencil=[-1，0]）#settingingtheherhsofthedifferentialequationDE.set_rhs（和dvetionTerm）expr（order，directtionName，time，stencil）········阿克斯4M. 卡拉姆，J.C.Sutherland和T.Saad/SoftwareX 12（2020）100541阿斯特丽德2i、ji−1，ji，j−1i−1，j− 1x.（十一）阿斯特丽德term2= （ DE.u （ time=n ， x=i ， y=j ） -DE.u（ time=n， x=i-1，y=j）-DE.u（ time=n，x=i， y=j-1） +DE.u （ time=n ， x=i-1， y=j-图1.一、使用PyModPDE 生成修改后的方程的步骤。term1=expr（order =1，directionName=回想一下，directionName只接受一个字符串值作为参数，这意味着不能使用此成员函数表示还有，时间给定术语的两个级别不能混合，这意味着表达式的离散表示中的所有元素必须在同一个timeIndex中进行计算。这意味着一个人不能在使用上述两种方法之一或两者的组合构建RHS表达式的项，用户可以通过调用成员函数set_rhs来设置微分方程的RHSset_rhs（expression）其中表达式是刚刚构造的项的线性组合，即，使用类似（u n+1−u n ）这种方法。set_rhs（term1+term2）i i−1第二种方法，以建立条款的RHS提供了更多的用户控制的离散化，并使一个构造任意复杂的有限差分格式与任意sten-numbers在时间和空间。它有效地模仿了我们手写有限差分公式的方式。这是通过使用<依赖变量名>（.）其具有以下签名（time，**kwargs）其中，时间是估计因变量的离散时间水平，例如，n+1，n，\ldots和kwargs是估计该表达式的空间点的索引，例如，x= i+1，y=j，\ldots.作为如何使用这种方法的一个例子，考虑方程中给出的离散化。（十）、然后，有限差分表示将转换为以下内容至此，设置完成，可以通过调用成员函数generate_modified_equation来生成修改后的公式generate_modified_equation（nterms）其中，nterms是正整数，其指示修改后的等式中要返回的项的函数generate_modified_equation（.）更改DifferentialEquation对象的状态，并允许使用辅助函数来访问各种形式的修改后的方程和放大因子。例如一可以调用display_modified_equation（）来显示修改后的公式的渲染LATEX表示。对于该示例DE. display_modified_equation（）将显示以下u与使用expr的第一种方法不同，这种方法可以用于表示交叉导数项。作为示例，考虑以下在时间t=tn处评估的二阶交叉导数项阿勒特 = −ax+2（x−ta）x2。（十二）此外，还可以通过调用DifferentialEquation对象的其他辅助函数来获取修改后的方程的其他可用形式例如，调用以下函数，乌鲁河nnn-un+un我们可以这样{\fn黑体\fs22\bord1\shad0\3aHBE\4aH00\fscx67\fscy66\2cHFFFFFF\3cH808080}{\fnarialblack\fs12\bord1\shad0\4aH00\fscx90\fscy110}{\fnarialx^{2}（DE. u（time=n，x=i，y=j）-DE. u（time=n，x=i-1，y=j））/DE. DxDE. latex_modified_equation（）=将返回一个LATEX字符串表示的修改后的方程，可以打印出来M. 卡拉姆，J.C.Sutherland和T.Saad/SoftwareX 12（2020）1005415= −==-我 =我（我 = −ci+1200xi−1.（十五）ut= − cux+tc uxx− cx+阿勒特阿克斯2X263）3号线e^{\alpha\Deltat}=1-a\frac{\Deltat}{\Deltax}\left（1-e^{-\ik_{x}\Deltax}\right）此外，用DE.symbolic_modified_equation（）实现了修正方程的符号形式.此外，一旦设置了微分方程的RHS，就可以调用成员函数generate_amp_factor（ ）来生成放大因子eαt。然后打电话给他。display_amp_factor（），用户可以显示放大因子的渲染LATEX形式。例如，放大因子的PDE离散方程。（9）被渲染和显示如下eαt1at阿克斯（1−e−ikx<$x），（13）这是该方案的预期放大因子上不同的注意，放大因子的帮助函数可以在调用DE之后使用。generate_modified_equation（.）而不需要先调用DE.generate_amp_factor（）。类似于修正后的方程，放大系数eαt的其他表示形式也可用。例如，可以调用函数DE. latex_amp_factor（）它将返回可以打印放大因子的LATEX字符串表示形式此 外 ， 放 大 因 子 的 符 号 形 式 也 是 可 用 的 ， 可 以 通 过 调 用DE.symbolic_amp_factor（）来访问。3. 说明性实例为了测试PyModPDE，我们研究了一维和二维时间相关PDE的修改方程的生成。3.1. 在一维平流上使用BTCS的示例作为上面提出的框架的一个例子，考虑平流方程ut=−cux;c>0，（ 14）使用时间上的向后欧拉离散化和空间上的中心离散化（BTCS）来求解。差分方程为un+1 − un（un+1−un+1）Δt图二. 在周期域上使用BTCS的高斯函数的平流，CFL=4。阿克斯0的情况。02，0的情况。01，并在CFL运行4.从第二在修改后的方程Eq. （18）我们知道，是扩散的，无论我们选择什么值作为振幅的减小，随着时间的推移，高斯函数的解是提前的。第三项在Eq。 公式（18）也是非零的，而与Δx和Δt的值无关，这表明该方案是色散的。 这种行为如图所示。由于与原始高斯函数相比，每个波数以不同的速度行进，因此2被称为为了使用PyModPDE获得该平流方程的修改方程，我们只需设置以下内容设un eat eikx.代入差分方程后，我们求解α1α=Δtlog2016年10月22日2xeixk+tce2ixk−tc（十六）最后，根据Eq。（4）我们有1′1′′BTC2a1=iα（0）= −c;a2=2i2α（0）=2;（十七）a=1′（0）= −c（nx2 +2c2）三点三！i3α六点三 ，修改后的方程的渲染乳胶表示为：显示为修正后的方程为二（二）263u=−c正如该计划所预期的那样。扩散项的系数为正，表明这是一个扩散格式，修正方程中的奇导数表明数值解中存在色散我们测试了这个数字-通过平流高斯函数u，u（x，t= 0）=e−x2/0。1，在-2到2的周期域上使用BTCS方案，3.2. 在2D平流上使用FTUS的示例作为多维示例，考虑二维中的平流u t= −bu x−cu y;（b，c）∈R+（19）fromsrc. pymodpde importDifferentiEquationfr e nfre nfrenfremport symbols#declare 为每个城市、空间和时间安排提供了一个系统，c，i，n=symbols（#constructatimede dentdifferental 等式DE=DiferentiEqation（depentVarName=“u”，independentVarsNames=[#我怎么能让我的孩子们advectionTerm=DE. expr（order =1，directionName=#设置DE中的DifferntialequationDE。set_rhs（-c*advectionTerm）#generatendisplayhemdifid 如果你不想让我知道，DE. generate_modified_equation（nterms =3）#displayhemdified equationDE.display_modified_equation（）u xxx+···。（十八）6M. 卡拉姆，J.C.Sutherland和T.Saad/SoftwareX 12（2020）100541∑JL我⏐⏐=⏐⏐y2Δt阿克斯埃什基X对于由Eq. （二十一）其中11表示在该区域中存在各向异性扩散，t= −cXJL⏐⏐在时间上使用前向欧拉法，在空间上使用UPWIND（FTUS）进行离散化，如下所示n+1−ni、jnnnnnni、jΔt= −bi、ji−1，jc阿克斯i、ji，j−1，（20）埃什基其中（i， j）是指二维网格上的位置（不要与本手稿中的虚单位和其他索引混淆在这种情况下，放大因子指数α为：α=1ln[1−b<$t（1−e−ikx<$x）−c<$t（1−e−iky <$y）]。 （二十一）在二维中，修改后的方程具有一般形式ut=a10ux+a 01uy+a 11uxy+a 20uxx+a 02uyy+ a12uxyy+a21uxxy+···∞=j，l≥0j+l≥1j+luajl（二十二）图3.第三章。 高斯函数在二维周期域上的平流。其中系数由下式给出：a=1nαyk=0（二十三）inj！l！∂k∂k⏐a10=1∂αk xk021= −b，a01=i∂αky=−c，（二十四）a11=1∂ α= −bct，修改后的公式的呈现乳胶形式显示为我 kxu102u102ux-y方向。数值扩散项也是非零的b2u系数由下式给出：+2（x−tb）x2。一个20021=2！I21=2！I2∂2α第二代∂2αk2B=2（x −C=2（y −（b）（同上）。（二十五）4. 影响PyModPDE为用户提供了一个高层次的抽象，以描述时间依赖的有限差分表示在这种情况下，修改后的方程的前几项是ut= −bux−cuy−bctuxy+b（x−tb）uxxC+2（y−tc）uyy+···（二十六）dent微分方程以及生成和返回修改后的公式的乳胶表示。该软件在犹他大学的数值方法课程中被用作辅助工具，用于介绍修正方程的概念。学生使用此工具帮助他们了解数字艺术效果和发生第三项在右边的方程。 （26）是一个各向异性项，导致在平流方向上的收缩和在运动的正交方向上的扩散，见图。3.第三章。使用PyModPDE，几行代码可以生成此数值方案的修改fromsrc. pymodpdeImportDifferentiEquationfrmpymp ort symbols#declare symbolsc，b，i，j，n=symbols（#constructatimededentPDEDE=DiferntialEquation（depentVarName=#我是如何创建这个文件的：advectionTerm1=-b*DE. expr（order=1，directionName=advectionTerm2=-c*DE. expr（order=1，directionName=#设置DE中的DifferntialequationDE。set_rhs（advectionTerm1+advectionTerm2）#generatendisplayhemdifid如果你不想让我知道，DE. generate_mdified_equation（nterms =2）#displayhemdified equationDE.display_modified_equation−k=0.M. 卡拉姆，J.C.Sutherland和T.Saad/SoftwareX 12（2020）1005417在对于域分辨率和时间步长的不同组合的某个数值解中。通过使用PyModPDE，学生将花费更多的时间来理解他们正在使用的数值方法，并减少调试时间。此外，学生将能够获得某些计划被认为是稳定的条件。我们不知道有类似于PyModPDE的公开软件。 我们相信，将PyModPDE提供给社区将有利于那些参与数值分析和计算物理的人-无论是在教育还是研究方面。5. 结论在这篇文章中，我们认为，有一个更容易的技术，推导出一个修改后的方程为一个给定的有限差分近似。首先由Chang[6]提出，并由作者独立得出，该技术是基于首先找到放大因子指数α，然后对波数求其导数。这些导数对应于修改后的方程的系数。我们还将该方法扩展到更高的维度，其中8M. 卡拉姆，J.C.Sutherland和T.Saad/SoftwareX 12（2020）100541在进行高维分析时，交叉导数常常被忽略。我们展示了多维修正方程如何导致交叉导数项，从而将各向异性效应引入平流。然后，我们提出的方法在Python软件PyModPDE中实现，该软件可以生成线性时间依赖偏微分方程的乳胶形式的修改后的方程。它是使用不同的时间积分器在多个空间维度的例子进行验证未来的工作将集中在支持高阶时间积分器，导致多个修改后的方程。竞合利益作者声明，他们没有已知的竞争性财务利益或个人关系，可能会影响本文报告的工作确认作者感谢犹他大学化学工程引用[1]Hirt CW.有限差分方程的启发式稳定性理论。 J计算Phys1968;2（4）：339-55.http ： //dx.doi.org/10.1016/0021-9991 （ 68 ） 90041-7 ，URLhttp://www.sciencedirect.com/science/article/pii/0021999168900417。[2]Klopfer G，McRae DS.分析有限差分格式的非线性修正方程法.第五届计算流体力学会议; 1981年。p. 1029号[3] Villatoro FR，Ramos JI.关于修正方程法第一章欧拉前向差分法的渐近分析 。 应 用 数 学 计 算 1999;103 （ 23 ） ： 111-39.http ：//dx.doi.org/10.1016/S0096-3003（98）10031-0，URLhttp://www.sciencedirect.com/science/article/pii/S0096300398100310。[4]莫莫尼亚特湖正则化长波方程非标准差分格式选择的修正方程法。 Abstr ApplAnal2014;2014：e754543.http://dx.doi.org/10.1155/2014/754543，URLhttps://www.hindawi.com/journals/aaa/2014/754543/abs/。[5]Warming RF，Hyett BJ. 修正后的方程更接近于有限差分法的稳定性和精度分析。JComputPhys1974;14（2）：159http://www.sciencedirect.com/science/article/pii/[6]张世春对修正的增温和海厄特方程技术的批判性分析。J Comput Phys1990;86（1）：107-26.[7] Ramshaw JD差分格式的数值粘性。Commun Numer Methods Eng 1994;10（ 11 ） ： 927-31.http://dx.doi.org/10.1002/cnm.1640101108 ， URLhttp://onlinelibrary.wiley.com/doi/10.1002/cnm. 1640101108/摘要。[8]Carpentier R，de La Bourdonnaye A，Larrouturou B.关于推导修 正 方程的线性分析数值方法。RAIROModel Math Anal Numer1997;31：459-70.