毕业设计（论文）-基于四次样条的矩阵微分方程数值解.doc

目录中文摘要：2abstract：21引言21.1几个重要名词的解释31.11函数矩阵31.12 函数矩阵的微分31.13样条函数42 常微分方程数值解法52.1 euler法：62.2法：82.3 线性多步法：92.4 外推法:133 基于四次矩阵样条的微分方程求解算法144 示例195 总结28致 谢28参考文献28基于四次样条的矩阵微分方程数值解中文摘要：矩阵微分方程经常出现在许多物理模型和工程技术模型中。利用四次矩阵样条构造形如，的一阶矩阵线性微分方程初值问题的数值解。给出实现算法和数值解的近似误差估计以及一些数值实例。并与利用三次矩阵样条构造一阶矩阵线性微分方程的数值解的结果进行了比较，比较结果表明用四次矩阵样条所构造的数值解的逼近效果要比用三次矩阵样条所构造的数值解的逼近效果要好。关键词：一阶线性矩阵微分方程，三次样条，四次样条，近似解abstract：matrix differential equations appear frequently in a wide variety of models in physics and engineering.this paper deals with the construction of approximate solution of first-order matrix linear differential equations giving by , using quartic cubic splines. an estimation of the approximation error, an algorithm for its implementation and some illustrative examples are included.compared with the results for the the same first-order matrix linear differential equations giving by , using cubic spline and quartic spline, the quartic spline is better than the cubic spline. keywords：first-order matrix linear differential equations, cubic spline, quartic spline.，approximation.1引言矩阵微分方程经常出现在许多物理模型和工程技术模型中1-4。除了许多数学模型中的问题经常写成矩阵的形式外，许多特殊的方法为了解决数量或向量问题也经常以矩阵形式出现。这方面的例子有许多，例如：用嵌入法求解具有边值条件的常微分方程的问题5；用打靶法求解具有边值条件的标量和向量分方程的问题6；用线性方法研究偏微分方程的数值积分问题7；用同伦方法求解非线性方程组的问题8。 向量化方法即将一个矩阵方程化为几个孤立的标量方程或向量方程有许多缺点9；。首先，用向量化方法会使方程中的量失去物理意义；其次，计算的成本会增加很多；最后，向量化方法没有利用一些符号语言能够处理矩阵表达式的优点。本文的组织结构如下：在第二部分中，将简要介绍常微分方程的一些常用解法，在第三部分中，我们给出了用四次矩阵样条来构造一阶矩阵微分方程近似解的方法，研究了该方法的逼近误差并编制了实现该方法的一个算法。最后在第四部分中给出了一些数值例子。在这篇论文里，我主要做了以下工作，在原有三次样条算法的基础上，利用四次矩阵样条构造一阶矩阵线性微分方程的数值解，并给出实现算法和数值解的近似误差估计以及一些数值实例。做这份论文的意义在于通过数值实例表明用四次样条构造一阶矩阵线性微分方程的数值解的逼近效果要比用三次矩阵样条所构造的数值解的逼近效果要好，或者说要得到相同的计算精度用四次样条比三次样条计算所需要的计算次数要少。1.1几个重要名词的解释1.11函数矩阵定义：如果矩阵a的每一个元素都是变量t的函数，则称a(t)= 为函数矩阵。1.12 函数矩阵的微分定义：设函数矩阵a(t)所有元素(t)都在点或某区间内可微，则称矩阵a(t)在点或某区间内是可微的。若a(t)可微，则其导数如下=()性质：设函数矩阵a(t)，b(t)都可微(1)若k为常数，则 (2)若a(t)，b(t)是同型矩阵，则(3)若a(t)是阶矩阵，b(t)是矩阵，则=特别地，如果a(t)或b(t)是常数矩阵a或b，就有，与数量函数微分法则不同的是，一般，即(4) 是t的可微函数1.13样条函数所谓样条即是在船体、汽车或航天器的设计中，模线设计员使用的弹性均匀的、窄的木条(或刚质条)。模线员在绘制模线时，用压铁压在样条的一批点上，强迫样条通过一组离散的型值点。当样条取得合适的形状以后，再沿着样条画出所需要的曲线，这将是一条光滑的曲线。样条函数最初就是来源于这样的样条曲线。用严格的数学定义表示就是下面形式定义：设在区间上给出一个划分：s(x)为实值函数，如果它满足条件(1)s(x)在子区间是n次多项式；(2)s(x)在区间上具有直到n-1阶连续微商即s(x)，则称y=s(x)为n次样条函数。当n=3时，称为三次样条函数；当n=4时，称为四次样条函数。2 常微分方程数值解法常微分方程的一般形式是 (2.1)这是一个n阶的常微分方程。如果给方程增加在x=x的一些条件，就变成了初值问题。另外的条件可以形成边值问题。微分方程最简单也是最重要的是一阶方程 (2.2)求满足(2.2)问题函数y(x)称为初值问题，我们要求解必须在区间x,b上确定。并且要求存在常数l对于x,b中的任何x与任何数，有 (2.3)这个条件称为条件。由常微分方程理论，如果满足条件，则在x,b上方程（2.2）存在唯一解。 虽然初值问题（2.2）对很大一类右端函数有解，但求出所需的解并非易事。实际上，除了极特殊情形外，人们不可能求出它的精确解，只是用各种近似方法得到满足一定精度的近似解。有一类近似解法只给出解在一些离散点上的值，称为数值方法。这一类方法的应用范围很广，特别适合用数字计算机计算。二阶方程与高阶方程可以化为一阶方程组。例如，对于二阶方程 (2.4)引入变量=z,则（2.4）可以写为 (2.5)因而高阶微分方程都等价于一阶微分方程组。一阶微分方程的常用解法有euler法，法，线性多步法，外推法等。下面将前面各解法简要回顾一下。2.1 euler法：这里，我们求一列点x,x,x上的近似解y,y.y,其中要求xxxx.xb (2.1.1)我们用等距节点，即h=x-x是常数，并称h是步长。在（2.2）中把x代入，则方程变成了 (2.1.2)为了求解，用差商y(x)-y(x)/h代替导数得 (2.1.3)用y(x)代替y代入式（2.1.3）,则有计算y的公式 (2.1.4)而y是（2.2）式中给出的初值，上述方法称为euler法，又称为euler折线法。由于用差商代替导数的方法有许多种，所以可以得到不同的euler法。用差商y(x)-y(x)/h代替,则由(2.2)得到 (2.1.9)同前，用用y代替y(x)，则有计算y的近似公式 (2.1.10)公式(2.1.10)称为隐式euler法。因为(2.1.10)式的右端含有未知量y，(2.1.10)式实质上是求y的一个方程。euler法(2.1.4)是计算y的一个直接公式，这类公式称为显式格式，公式(2.1.10)与公式(2.1.4)一样也是一个一阶公式，但是隐式公式。为了提高euler法计算格式的精度，我们采用中心差商代替导数的方法y(x)-y(x)/2h代替，则由等式(2.1.2)，且用y代替y(x),得计算y的近似公式 y=y+2h f(x,y) (2.1.11)公式(2.1.11)是两步euler格式。在式子(2.1.4)与式子(2.1.10)中，计算y只用到前一步y的信息，所以它们都是一步法。公式(2.1.11)计算y时，既要用到前一步的信息y,又要用到y，所以称为两步法。在(2.1.11)式的计算中，我们已知初值y,可以用(2.1.4)或(2.1.10)计算，当然也可用后边得到的高阶公式求得。 由数值微分公式 f(x)=- (2.1.12)所以euler两步格式(2.1.11)局部误差为 所以euler两步格式(2.1.11)是二阶方法。其中方法(2.1.11)局部误差为h(x,h)=y(x)+2hf(x,y(x)-y(x+h)。也可从另一个角度来解微分方程（2.2）.对方程在小区间（x，x）上积分，得y(x)-y(x)= (2.1.14)为计算（2.1.14）右边的积分，我们使用梯形公式，得所以y(x)y(x)+ f(x,y(x)+ f(x,y(x) (2.1.15)同前，用y(x),的近似值y代入（2.1.15）式，则由计算y的公式，称为梯形格式 y= y+ f(x,y(x)+ f(x,y(x) (2.1.16)注意梯形公式的误差估计- f(x,y(x)+ f(x,y(x)=-() (2.1.17)则梯形格式（2.1.16）的局部误差为h(x,h)= () (2.1.18)所以梯形格式也是二阶方法。事实上，它是euler显式格式（2.1.4）与euler隐式格式（2.1.4）的算术平均，但方法的阶却提高了。 由前面知道，euler法是一个显式格式，它好计算，是一步法，且计算量少，但只是一个一阶方法，计算精度比较差。而上面得到的梯形方法也是一步的方法，而是二阶方法，但是一种隐式格式，需要解关于y的方程（通常的方法是迭代法），计算量大。 综合上述两种方法的优点，就得到了一种新的算法，就是下面的法。2.2法：法是按关系式y= y+h (2.2.12)逼近y(x)，其中k满足 k=f(x+ y+h) (2.2.13)在式子(2.2.12)、(2.2.13)中，,都是常数为0以便使k= f(x,y).,可以用待定系数法确定，以得到比较高阶的公式。由式子(2.2.13)给出k的公式(2.2.12)是隐式公式，显然法写为 y= y+hk=f(x+ y+h) (2.2.14)它可以递归的给出k的值而不需要解非线性方程组法的实质是用h去逼近积分 (2.2.15) 用此思想，我们把式(2.2.12)用下述形式表示。设x。以x代替x，以y代替y，其中y是一任意实数。引进函数 (2.2.16)则(2.2.12)可以写为简洁的形式 (2.2.17)上式形式上看与显式euler法类似，不同的是这儿用函数代替了.这时，法(2.2.12)、(2.2.13)与(2.2.14)在点x的局部误差为 (2.2.18) 其中y(x)表示初值问题的精确解2.3 线性多步法：前边的方法是计算常微分方程初值问题的一类重要的算法。但是，这类算法在每进行一步都需要先计算几个点上的斜率值（即k, k,）,计算量比较大，并且，在计算y之前，我们已经得到了y，y，的值，而方法只用到了y，而没有用到前几步计算得到的信息。现在，考虑点,0jn,其中n=(b-)/h.假定近似函数值y，y，.,y,y已知（实际上，开始值可由方法求得）。解的导数可用,k+1-ijk来估计。线性多步法的一般公式是 (2.3.1)其中系数，,是固定的数。0时y的值也在右端项中出现。故在这种情况下决定必须解由（2.3.1）得到的非线性方程，这时（2.3.1）式是线性多步法的隐式公式，为强调右端对的依赖性将它重新写成 (2.3.2)有前面引进的条件的假定下，对于充分小的值h，可以证明，方程（2.3.2）可用简单迭代法求解。设表示方程（2.3.2）的右边，对任意值与当步长h取得足够小，使1时，（2.3.2）右端的算子是压缩的，所以可以用非线性的方法解得。例如可用简单迭代与迭代等方法。当=0时，得到线性多步法的显式公式假定函数z(x)是可微的，并且引进差分算子lz= (2.3.4)其中取，任意lz=0等价与没有截断误差的函数z(x)的线性多步法的公式。由线性平移和自变量的尺度变换，特别可取h=1,k=和x=0,这样算子化为lz= (2.3.4)现在，线性多步法的阶由使截断误差为零的最高次多项式的次数来定义。 定义7.3.1 如果p是使 lx=0, =0,1,2,p (2.3.5)的最大正整数，则称线性多步法(2.3.1)是p阶的。 为了方便，记0=1，这时方程(2.3.5)可重新写为=0 (2.3.6)其中(-1)：=max0, -1.因=1,我们得到关于2i+1个未知数的p+1个方程。显然，当p2i时，(2.3.6)式无解, 当p=2i时，(2.3.6)式有唯一解。 定义2.3.2 设y是方程 ， y(x)=y的精确解，方法(2.3.1)的局部误差由 h=h(x,h)=ly (2.3.7)定义，其中l由(2.3.4)式给出。 假定y是解析的，可以容易的算出h的值。对于任何整数，利用taylor公式y(x+h)= y(x)+ (2.3.8)对于p阶方法，得=于是局部误差的阶和方法的阶相同。 对于不同的p和值，解方程(2.3.6)后，就得到不同标准的方法。 对p=2和=1时，方程(2.3.6)变成 对于=0，1+=0 对于=1，1-=0 对于=2，1-2=0由此得到=-1，=1/2，=1/2。这时对应的公式是=y+ (2.3.9)这就是euler法中的隐式格式，有时也称为修正的euler公式。重复应用公式(7.3.3),多步公式的局部误差和舍入误差都会有积累，可以证明，当且仅当方程 p(z)=0 (2.3.10)在单位圆外无零点，在单位圆周上零点重数为1时，对于小的h值和小的舍入误差，则舍入误差和初值函数的误差函数的误差传播影响也是小的。具有这种性质的方法称为稳定的方法。多步法的计算问题： 用步法计算时，须知道初值，其中y(x)=y是给定的初值，其余是附加初值。计算附加初值有几种方法。一种是单步法，比如euler法和法及其他单步法。为了保持多步法的精度，计算附加初值时要将x,x之间节点加密或采用和多步法有同样阶的法。 多步法的第二个问题是如何选择阶p(或者步数k)。从收敛阶的观点，自然希望把p取大一些。但是高阶收敛方法要求解的光滑型也很高，否则达不到高精度的密度。并且高阶多步法的绝对稳定域也小，所以p的大小要考虑到解的光滑型和稳定性以及总的工作量。 多步法的第三个问题是步长h的选取。理论上似乎按照误差估计式选定h是合理的，但那种估计往往偏大，因此选定的h可能过小，既不必要也不经济。实际用的步长h不是一次取定，而是根据精度要求，又粗到细逐渐调整（选步长），当h达到要求后就以此为步长计算。在计算中还可以改变步长，但计算过程变复杂了，这就是多步法的缺点。与此相反，单步法适合用变步长计算2.4 外推法:在求解初值问题 (2.4.1)的时候，可以写为 (2.4.2)而(2.4.2)式右端第二项的积分可以用各种各样的数值积分求解。这当然也可以用外推法来求（2.4.1）的解。让我们回忆一下richardson外推法。设t是某个被逼进的对象，t(h)是依赖于步长h的t的一个近似式 tt(h) (2.4.3)如果t与t(h)之间能表示为渐近式 (2.4.4)而0,且是与h无关的常数。这时richardson外推法计算格式为 (2.4.5)所以使用外推法的关键，是看存在什么样的渐近展开式(2.4.4)。对于euler法，我们知道局部截断误差可写为 (2.4.6)因此，经过几次计算得到后，鉴于全程截断函数误差比局部截断误差的阶数低一次，可设想 (2.4.7)对于梯形法，因为它是由梯形求积公式推出来的，而梯形求积公式的截断误差具有展开式所以由前面的知识可知梯形法的局部截断误差为 (2.4.8)如前，经过几次计算得到后，可设想 (2.4.9)同样，对于由中心商得到的euler中心法公式 y=y+2h f(x,y) (2.4.10)(2.4.8)也成立。根据(2.4.5)，令，外推法计算公式为 (2.4.11)其中，对(2.4.7)，=k, =1,对(2.4.9)，=2k, =2,。注意是步长为h时求得的的近似值（例如为在步长为h时的近似值，这时就取=）。与前面不同的是=。3 基于四次矩阵样条的微分方程求解算法本文的工作是给出如(1)所示的一阶线性矩阵微分方程的一种数值求解方法 y(x) = a(x)y(x) + b(x), a x b, (1)其中,yc,a: c,b: c满足a,bc(a,b),它保证了（1）只有一个解，并且这个解是连续的、可微的 10，p99。在参考文献11中，有许多不同的工程应用问题的数学模型都是以（1）的形式出现的。一些非线性的方程（如黎卡提方程12-14）通过线性化后也可以化为（1）的形式。利用固定步长的线性多步法来解决诸如（1）类型的方程的近似数值方法在参考文献15中有详细介绍，但是这个误差的范围依赖于步长h，且是h的若干次指数次幂，所以在实际操作中应使h尽可能小，并且在这些方法中为了得到一个连续的解，还必须要采用一些插值技巧。另外一些方法则是从magnus or fer方法发展起来的，利用这些方法求解问题（1）时，要计算一些矩阵的的幂，这就增加了计算成本。 文献17给出了用三次样条构造标量形式的常微分方程的数值解的方法，所得近似解具有以下优点：1、近似解在区间a , b上是连续的；2、近似解易于求值计算；3、近似解所达到的误差阶是。最近，文献18也用了样条方法来构造另外一些标量型的微分方程的近似数值解。参考文献19研究了利用基于hermite插值技巧的隐式样条方法来构造向量型的微分方程的数值解的问题。在本文中，我们给出了利用四次矩阵样条来构造一阶矩阵微分方程（1）的近似数值解的方法。并与利用三次矩阵样条构造一阶矩阵线性微分方程的数值解的结果进行了比较，比较结果表明用四次矩阵样条所构造的数值解的逼近效果要比用三次矩阵样条所构造的数值解的逼近效果要好。所得方法仍然具有参考文献17中方法所具有的优点。 在这次工作中我们以c表示p阶复矩阵，如果ac我们定义矩阵a的2次范数为=，在这里z是一个属于的复向量,=(zz)是z的欧几里得范数，根据参考文献20，p.56，我们可得到满足如下形式的不等式 (2)我们以p(x)表示以实变量x为自变量的n次矩阵多项式。设g是定义在区间a,b上的一个矩阵函数，即g: c,如果g在区间上是k次可微分的，且它的第k阶导函数在上是连续的，则称g是一个k类。记为gc()。取r和, 是上的一个划分，给定一个大于等于0的整数m,我们定义m次矩阵样条函数构成的集合如下： 如果m=3，则该矩阵样条称为三次矩阵样条。如果m=4，则该矩阵样条称为四次矩阵样条。让我们给出区间上的一个划分,， (3)其中h=(b-a)/n, n是一个确定的整数。 在问题（1）中我们构造一个固定步长的四次矩阵样条作为近似解，对于第一个间隔,我们把样条函数定义为如下形式， (4)其中矩阵 c是一个待定参数，（4）中的s(x)满足s(a)=y(a) , (a)=(a)=a(a)y(a)+b(a),为了得到样条函数我们必须知道(a), (a)和，利用a、b经过简单计算可得到 (5)利用(5)可以得到(a)。利用上式继续对x求导可得(x)= a(x)+a(x)+ y(x)+a(x)(x)+ a(x)+ b(x)+ (*)从而可以得到 我们还只剩下唯一的没有解出，一旦解出那么样条函数就可确定。因而，我们可以利用样条函数在这一点上满足微分方程(5)即 (6)其中在(6)中我们得到一个只包含一个未知矩阵参数的矩阵方程即(7) 假设矩阵方程(7)有解，那么利用这个方法就得到了一个定义在上的样条函数。在区间a+h,a+2h上，四次矩阵样条的函数形式为s(x)= s+ + +, (8)所以定义在区间上的s(x)是a,a+2h上的2阶可导函数并且样条函数s(x)只有一个未知系数 c，容易验证样条函数(8)在点上满足微分方程（1），同样，样条函数(8)在点x=a+2h处也满足微分方程（1），从而可以确定即 (9)化简(9)得到一个只包含一个未知矩阵系数的矩阵方程即 (10)如果矩阵方程(10)只有一个解，这样定义在上的样条函数就可以完全确定。重复这一过程，假设我们已经得到定义在上的四次矩阵样条函数，那么我们将把定义在区间上的样条函数定义为s(x)=(x)+, (11)其(x)=s (12)至此，我们就得到了样条函数,容易验证它在点处满足微分方程（1），又样条函数在点处满足微分方程（1）从而可以得到即经过化简可以得到矩阵微分方程 (13)注意由于矩阵(13)是可逆的所以方程(13)的可解性就得到了保证，其中k=0,1n-1.如果用m=max;ax0且,m由(14)式给出，那么形如(11)所示的四次样条函数s(x)在每一个小区间 (k=0,1.n-1)上就被完全定义下来并且只能是此种形式；(2) 如果a,b c(),那么=o(h),x下面的算法就是利用定义在区间上四次矩阵样条得到的问题(1)的近似解。算法：step1：利用等式(14)和(5)分别确定常数m和(a)，取nm(b-a)/4， 以及如(3)所示给出的上的一个划分；step2：对于k=0解矩阵等式(7)并且计算出如(4)给出的s(x)；step3：对于k=1,.n-1解矩阵等式(13)并且计算出如(11)给出的s(x)。4 示例在这一部分，我们主要是测试该算法在实际计算过程中的精确程度及算法的稳定性并于三次样条计算得到的结果作比较。例4.1 考虑问题 y(x) = y(x)， 0x1, , y(x) c, (16)这一问题的已知解是，所以我们可以计算出近似误差。对(16)求导可得(x)=(+a(x)y(x)， y(0)=a(0)y(0)= ，(0)= ，并且因为，所以我们取m=3,对于n=10,它满足nm(b-a)/4.所以我们得到h=(b-a)/n=0.1.利用mathematicav4.0可以得到结果如表1所示。用三次样条逼近得到的结果如表2 所示 。表1近似解近似误差表2近似解近似误差11.536611.557916.36816.404723.306324.7135比较三次样条和四次样条的计算误差，可以看出四次样条的计算误差是三次样条的十分之一，计算精度比三次样条提高许多。对于相同的数值逼近问题，如果取相同的步长，四次样条的精度比三次样条高一个数量级，也就意味着如果要取得相同的计算精度，四次样条的计算次数只是三次的一半，所以利用四次样条可以减少计算量。 例4.2对于矩阵问题 其中 并且该问题的精确解是 因为所以我们取m=5，由此得到下式(x)= a(x)+a(x)+ y(x)+a(x)(x)+ a(x)+ b(x)+ ，并且可得 ，对于n=10,它满足，所以。利用前面的算法同样可以得到用三次样条和四次样条逼近的对比结果，其中四次样条逼近的结果如表3所示，三次样条的逼近结果如表4所示。 表3近似解近似误差0，0.1表4近似解近似误差0，0.1比较上面两表中数据，同样可以看出用四次样条比用三次样条的计算精度高，四次样条的逼近效果更好。5 总结在数值逼近问题中，用四次样条的逼近效果比三次样条要好。四次样条的计算误差只是三次样条的十分之一，计算精度比三次样条提高许多。对于相同的数值逼近问题，如果取相同的步长，四次样条的计算精度比三次样条高一个数量级；如果要取得相同的计算精度，四次样条的计算次数只是三次的一半，所以利用四次样条可以减少了计算量。 致 谢非常感谢郭清伟老师在我论文写作过程中的鼎力帮助，感谢陆琴、邢芳等同学对于我论文提出了宝贵的意见，和所有合肥工业大学的老师对我的帮助。参考文献i. l.d. faddeyev, the inverse problem in the quantum theory of scattering, j. math. physics 4 (i), 72-104,(1963).2. s. barnett, matrices n control theory, van nostrad, reinhold, (1971).3. w.t. reid, riccat differential equations, academic press, (1972).4. r.a. labudde, classical mechanics of molecular collisions, university of wisconsin theretica chemistry inst. report wis-tc1-414, (1973).5. m. scott, inariant mbeddizg and is applications to ordinary differential equations, addison-wesley,(1973).6. p. marzulli, global error estimates for the standard parallel shooting method, j. computat. appl. mahs.34, 233-241, (1991).7. k. raktorys, the method of dicre$zation in time and partial differential equations d. reidel pub. co.,dordrecht, (1982).8. p.t. boggs, the solution of nonlinear systems of equations by a-stable integration techniques, siam j.numer. anal. 8 (4), 767-785, (1971).9. a. graham, kvanecker producs and matrix calculu with applications, john wiley, new york, (1981).10. t.m. flett, differential analysis, cambridge university press, (1980).11. u.m. ascher, r.m.m. mattheij and r.d. russell, numerical solutions of boundary value proble for ordinary differential eations, penticc hall, new jersey, (1988).12. l. jsdar and e. ponsoda, non-autonomous riceati-type matrix differential equations: existence interval,construction of continuous numerical solutions and error bounds, ima j. numer. anal. 15 (1), 61-74,(1995).13. l. jddar, j.c. cortds and j. l. morsra, construction mad computation of variable coefficient sylvester differential problems, computers maths. applic. 32 (8), 41-50, (1996).14. l. jddax and j.c. corts, rational matrix approximation with a priori erro