基于四次样条的矩阵微分方程数值解

目录

中文摘要： ............................................... 2 Abstract： ............................................... 2 1引言 ................................................... 3

1.1几个重要名词的解释 ................................................................................................. 3

1.11函数矩阵 .......................................................................................................... 3 1.12 函数矩阵的微分 .............................................................................................. 4 1.13样条函数 .......................................................................................................... 4

2 常微分方程数值解法 ..................................... 5

2.1 Euler法： ................................................................................................................... 6 2.2Rung-Kutta法： ...................................................................................................... 8 2.3 线性多步法： ............................................................................................................ 9 2.4 外推法: ..................................................................................................................... 13

3 基于四次矩阵样条的微分方程求解算法 .................... 14 4 示例 .................................................. 19 5 总结 .................................................. 28 致 谢 ................................................... 28 参考文献 ................................................ 28

第1页

基于四次样条的矩阵微分方程数值解

基于四次样条的矩阵微分方程数值解

中文摘要：矩阵微分方程经常出现在许多物理模型和工程技术模型中。利用四

次矩阵样条构造形如y'?A(x)y?B(x)，

y(0)?y0,x?[a,b];A(x),B(x)?C4[a,b]

的一阶矩阵线性微分方程初值问题的数值解。给出实现算法和数值解的近似误差估计以及一些数值实例。并与利用三次矩阵样条构造一阶矩阵线性微分方程的数值解的结果进行了比较，比较结果表明用四次矩阵样条所构造的数值解的逼近效果要比用三次矩阵样条所构造的数值解的逼近效果要好。

关键词：一阶线性矩阵微分方程，三次样条，四次样条，近似解

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

y'?A(x)y?B(x),y(0)?y0,x?[a,b];A(x),B(x)?C4[a,b]

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

y'?A(x)y?B(x),y(0)?y0,x?[a,b];A(x),B(x)?C4[a,b]

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.

第2页

1引言

矩阵微分方程经常出现在许多物理模型和工程技术模型中[1-4]。除了许多数学模型中的问题经常写成矩阵的形式外，许多特殊的方法为了解决数量或向量问题也经常以矩阵形式出现。这方面的例子有许多，例如：用嵌入法求解具有边值条件的常微分方程的问题[5]；用打靶法求解具有边值条件的标量和向量分方程的问题[6]；用线性方法研究偏微分方程的数值积分问题[7]；用同伦方法求解非线性方程组的问题[8]。

向量化方法即将一个矩阵方程化为几个孤立的标量方程或向量方程有许多缺点[9]；。首先，用向量化方法会使方程中的量失去物理意义；其次，计算的成本会增加很多；最后，向量化方法没有利用一些符号语言能够处理矩阵表达式的优点。

本文的组织结构如下：在第二部分中，将简要介绍常微分方程的一些常用解法，在第三部分中，我们给出了用四次矩阵样条来构造一阶矩阵微分方程近似解的方法，研究了该方法的逼近误差并编制了实现该方法的一个算法。最后在第四部分中给出了一些数值例子。

在这篇论文里，我主要做了以下工作，在原有三次样条算法的基础上，利用四次矩阵样条构造一阶矩阵线性微分方程的数值解，并给出实现算法和数值解的近似误差估计以及一些数值实例。做这份论文的意义在于通过数值实例表明用四次样条构造一阶矩阵线性微分方程的数值解的逼近效果要比用三次矩阵样条所构造的数值解的逼近效果要好，或者说要得到相同的计算精度用四次样条比三次样条计算所需要的计算次数要少。 1.1几个重要名词的解释 1.11函数矩阵

定义：如果矩阵A的每一个元素aij都是变量t的函数，则称

?a11(t)...a1n(t)???A(t)= ?...................?为函数矩阵。

??a21(t)...ann(t)??第3页

基于四次样条的矩阵微分方程数值解

1.12 函数矩阵的微分

定义：设函数矩阵A(t)所有元素aij(t)都在t0点或某区间内可微，则称矩阵A(t)在t0点或某区间内是可微的。若A(t)可微，则其导数如下

ddA(t)=(aij(t)) dtdt性质：设函数矩阵A(t)，B(t)都可微

(1)若k为常数，则 [kA(t)]'=kA'(t)

(2)若A(t)，B(t)是同型矩阵，则[A(t)+B(t)]'=A'(t)+B'(t) (3)若A(t)是m?n阶矩阵，B(t)是n?s矩阵，则

[A(t)B(t)]'=A'(t)B(t)+A(t)B'(t)

特别地，如果A(t)或B(t)是常数矩阵A或B，就有

[AB(t)]'=AB'(t)，[A(t)B]'=A'(t)B，与数量函数微分法则不同的是，一般 [A(t)B(t)]'?A'(t)B(t)+B'(t)A(t)，即[A2(t)]'?2A(t)A'(t)

(4)[a(t)A(t)]'=a'(t)A(t)+a(t)A'(t) a(t)是t的可微函数 1.13样条函数

所谓样条即是在船体、汽车或航天器的设计中，模线设计员使用的弹性均匀的、窄的木条(或刚质条)。模线员在绘制模线时，用压铁压在样条的一批点上，强迫样条通过一组离散的型值点。当样条取得合适的形状以后，再沿着样条画出所需要的曲线，这将是一条光滑的曲线。样条函数最初就是来源于这样的样条曲线。用严格的数学定义表示就是下面形式 定义：设在区间[a,b]上给出一个划分：

?:={- ?=x0 < x1 <...

(1)s(x)在子区间[xi,xi?1](i=0,1...,n)是n次多项式； (2)s(x)在区间[a,b]上具有直到n-1阶连续微商

第4页

即s(x)?Cn?1(a,b)，则称y=s(x)为n次样条函数。 当n=3时，称为三次样条函数； 当n=4时，称为四次样条函数。

2 常微分方程数值解法

常微分方程的一般形式是

y(n)=f (x,y,y',y'',.......y(n?1)) (2.1) 这是一个n阶的常微分方程。如果给方程增加在x=x0的一些条件，就变成了初值问题。另外的条件可以形成边值问题。

微分方程最简单也是最重要的是一阶方程

?y'?f(x,y) (2.2) ?y(x)?y00?求满足(2.2)问题函数y(x)称为初值问题，我们要求解必须在区间[x0,b]上确定。并且要求存在常数L对于[x0,b]中的任何x与任何数y1，y2有

|f(x,y )-f(x,y )|?dL|y -y | (2.3)

这个条件称为Lipschitz条件。由常微分方程理论，如果f(x,y)满足Lipschitz条件，则在[x0,b]上方程（2.2）存在唯一解。

虽然初值问题（2.2）对很大一类右端函数有解，但求出所需的解并非易事。实际上，除了极特殊情形外，人们不可能求出它的精确解，只是用各种近似方法得到满足一定精度的近似解。有一类近似解法只给出解在一些离散点上的值，称为数值方法。这一类方法的应用范围很广，特别适合用数字计算机计算。

二阶方程与高阶方程可以化为一阶方程组。例如，对于二阶方程

y''=f(x,y , (2.4) 引入变量y'=z,则（2.4）可以写为

y'?z? ? (2.5)

z'?f(x,y,z)?第5页

基于四次样条的矩阵微分方程数值解

因而高阶微分方程都等价于一阶微分方程组。一阶微分方程的常用解法有Euler法，Rung-Kutta法，线性多步法，外推法等。 下面将前面各解法简要回顾一下。 2.1 Euler法：

这里，我们求一列点x1,x2,…,xn上的

y1,y2…….yn,

其中要求x0

y'(xk )=f(xk ,y(xk )) (2.1.2)

为了求解，用差商[y(xk?1)-y(xk)]/h代替导数y'(x )得

y(x k+1) y(xk )+h f(xk ,y(xk )) (2.1.3)

用y(x1)代替y1代入式（2.1.3）,则有计算yk?1的公式

yk+1 = yk +hf(xk ,y k)k=0,1,2... (2.1.4) 而y0是（2.2）式中给出的初值，上述方法称为Euler法，又称为Euler折线法。

由于用差商代替导数的方法有许多种，所以可以得到不同的Euler法。 用差商[y(xk?1)-y(xk?1)]/h代替y'(xk?1) ,则由(2.2)得到

y(xk+1 ) ?y(xk )+hf(x k+1,y(x k+1)) (2.1.9) 同前，用用yl代替y(xl)，则有计算yk?1的近似公式

y k+1=yk +hf(x k+1 ,y(x k+1 )) (2.1.10) 公式(2.1.10)称为隐式Euler法。因为(2.1.10)式的右端含有未知量yk?1，(2.1.10)式实质上是求yk?1的一个方程。Euler法(2.1.4)是计算yk?1的一个直接公式，这类公式称为显式格式，公式(2.1.10)与公式(2.1.4)一样也是一个一阶公式，但是隐式公式。

第6页

为了提高Euler法计算格式的精度，我们采用中心差商代替导数的方法[y(xk?1)-y(xk?1)]/2h代替y'(x )，则由等式(2.1.2)，且用yl代替y(xl),得计算yk?1的近似公式 yk?1=yk?1+2h f(xk,yk) (2.1.11) 公式(2.1.11)是两步Euler格式。在式子(2.1.4)与式子(2.1.10)中，计算yk?1只用到前一步yk的信息，所以它们都是一步法。公式(2.1.11)计算yk?1时，既要用到前一步的信息yk,又要用到yk?1，所以称为两步法。在(2.1.11)式的计算中，我们已知初值y0,y1可以用(2.1.4)或(2.1.10)计算，当然也可用后边得到的高阶公式求得y1。 由数值微分公式

y(xk?1 )-y(xk-1 )h2f??(?k) (2.1.12) f(xk)=-2h6所以Euler两步格式(2.1.11)局部误差为

h3H(xk ,h)= ?f??(?k)

3所以Euler两步格式(2.1.11)是二阶方法。其中方法(2.1.11)局部误差为

H(xk,h)=y(xk?1)+2hf(xk,y(xk))-y(xk+h)。

也可从另一个角度来解微分方程（2.2）.对方程y'=f(x,y) 在小区间（xk，xk?1）上积分，得

y(xk?1)-y(xk)=?xk?1xkf(x,y)dx (2.1.14)

为计算（2.1.14）右边的积分，我们使用梯形公式，得

?xk?1xkhf(x,y?x?)dx? [ f(xk ,y(xk ))+ f(xk+1 ,y(xk+1 ))]

2h[ f(xk,y(xk))+ f(xk?1,y(xk?1))] (2.1.15) 2所以

y(xk?1)?y(xk)+

同前，用y(xl),的近似值yl代入（2.1.15）式，则由计算yk?1的公式，称为

第7页

基于四次样条的矩阵微分方程数值解

梯形格式 yk?1= yk+

h[ f(xk,y(xk))+ f(xk?1,y(xk?1))] (2.1.16) 2注意梯形公式的误差估计

?xk?1xkhh3f''(?k) (2.1.17)f(x,y)dx-[ f(xk,y(xk))+ f(xk?1,y(xk?1))]=-212h3则梯形格式（2.1.16）的局部误差为H(xk,h)= f''(?k) (2.1.18)

12所以梯形格式也是二阶方法。事实上，它是Euler显式格式（2.1.4）与Euler隐式格式（2.1.4）的算术平均，但方法的阶却提高了。

由前面知道，Euler法是一个显式格式，它好计算，是一步法，且计算量少，但只是一个一阶方法，计算精度比较差。而上面得到的梯形方法也是一步的方法，而是二阶方法，但是一种隐式格式，需要解关于yk?1的方程（通常的方法是迭代法），计算量大。

综合上述两种方法的优点，就得到了一种新的算法，就是下面的Rung-Kutta法。

2.2Rung-Kutta法：

Rung-Kutta法是按关系式yk?1= yk+h?aiKi (2.2.12)

i?1m逼近y(xk?1)，其中Ki满足

Ki=f(xk+?ih, yk+h??ijKj) i=0,1,...,m (2.2.13)

i?1m在式子(2.2.12)、(2.2.13)中，?i,?i,?ij都是常数?1j为0以便使K1= f(xk,yk).?i,?i,?ij可以用待定系数法确定，以得到比较高阶的公式。由式子(2.2.13)给出Ki的公式(2.2.12)是隐式公式，显然Rung-Kutta法写为

(j? i, ?ij =0) yk?1= yk+h?aiKi

i?1m第8页

Ki=f(xk+?ih, yk+h??ijKj) i=0,1,...,m (2.2.14)

i?1m它可以递归的给出Ki的值而不需要解非线性方程组Rung-Kutta法的实质是用h?aiKi去逼近积分?i?1mxk?1xkf(x,y?x?)dx (2.2.15)

用此思想，我们把式(2.2.12)用下述形式表示。设x?[x0,b?h]。以x代替xk，以y代替yk，其中y是一任意实数。引进函数

g(x,y,h)= ?aiKi (2.2.16)

i?1m则(2.2.12)可以写为简洁的形式

yk+1 = yk +hg(x,y,h) (2.2.17) 上式形式上看与显式Euler法类似，不同的是这儿用函数g(x,y,h)代替了 f(x,y).这时，Rung-Kutta法(2.2.12)、(2.2.13)与(2.2.14)在点x的局部误差为

H(x,h)=y(x)+ hg(x,y(x),h)-y(x+h) (2.2.18)

其中y(x)表示初值问题的精确解

2.3 线性多步法：

前边的Rung-Kutta方法是计算常微分方程初值问题的一类重要的算法。但是，这类算法在每进行一步都需要先计算几个点上的斜率值（即K1, K2,…）,计算量比较大，并且，在计算yk?1之前，我们已经得到了yk，yk?1，…的值，而Rung-Kutta方法只用到了yk，而没有用到前几步计算得到的信息。

现在，考虑点xj?x0?jh,0≤j≤N,其中N=(b-x0)/h.假定近似函数值yk?1?i，yk?2?i，….,yk?1,yk已知（实际上，开始值y1,y2...yi?1可由Rung-Kutta方法求得）。解的导数可用fi?f(xj,yj),k+1-i≤j≤k来估计。线性多步法的一般公式是

第9页

基于四次样条的矩阵微分方程数值解

yk?1???ajyk?1?j?h??jfk?1?j (2.3.1)

j?1j?0ii其中系数aj和?j，1?j?i,是固定的数。?0?0时yk?1的值也在右端项

fk?1?f(xk?1,yk?1,)中出现。故在这种情况下决定yk?1必须解由（2.3.1）得到的非线性方程，这时（2.3.1）式是线性多步法的隐式公式，为强调右端对yk?1,的依赖性将它重新写成

yk?1???ajyk?1?j?h??jfk?1?j?h?0f(xk?1,yk?1) (2.3.2)

j?1j?0ii有前面引进的Lipschitz条件的假定下，对于充分小的值h，可以证明，方程（2.3.2）可用简单迭代法求解。设g(yk?1)表示方程（2.3.2）的右边，对任意值yk?1,与yk?1

?|g(yk?1)-g(yk?1)|?h?0|f(xk?1,yk?1)?f(xk?1,yk?1)|?h?0L|yk?1?yk?1|

当步长h取得足够小，使h?0L<1时，（2.3.2）右端的算子是压缩的，所以可以用非线性的方法解得。例如可用简单迭代与Ssidel迭代等方法。

当?0=0时，得到线性多步法的显式公式

ii???yk?1???ajyk?1?j?h??jfk?1?j

j?1j?0假定函数z(x)是可微的，并且引进差分算子

iLkz=?[ajz(x0?(k?1?j)h)?h?jz?(x0?(k?1?j)h)] (2.3.4)

j?0其中取?0 =1，任意Lkz=0等价与没有截断误差的函数z(x)的线性多步法的公式。由线性平移和自变量的尺度变换，特别可取h=1,k=i和x1=0,这样算子化为

第10页

Lz=?[ajz(i?j)??jz?(i?j)] (2.3.4’)

?ij?0现在，线性多步法的阶由使截断误差为零的最高次多项式的次数来定义。 定义7.3.1 如果p是使 L?xl=0, l=0,1,2,…,p (2.3.5) 的最大正整数，则称线性多步法(2.3.1)是p阶的。 为了方便，记00=1，这时方程(2.3.5)可重新写为

i?[aj?0j(i?j)l??jl(i?j)(l?1)0]=0 l=0,1,2,...,p (2.3.6)

其中(l-1)0：=max{0, l-1}.因?0=1,我们得到关于2i+1个未知数的P+1个方程。显然，当p<2i时，(2.3.6)式有无穷多解，当p>2i时，(2.3.6)式无解, 当p=2i时，(2.3.6)式有唯一解。

定义2.3.2 设y是方程

y'(x)=f(x,y(x))， y(x0)=y0

的精确解，方法(2.3.1)的局部误差由 Hk=H(xk,h)=Lky (2.3.7) 定义，其中Lk由(2.3.4)式给出。

假定y是解析的，可以容易的算出Hk的值。对于任何整数?，利用Taylor公式

y?(xk?1)yp(xk?1)(?h)+…+y(xk?1+?h)= y(xk?1)+(?h)p+… (2.3.8) 1!p!对于p阶方法，得

yp?1(xk?1)p?1?p?1Lk y = hLt?O(hp?2)

(p?1)!=

yp?1(xk?1)p?1ih?[aj(i?j)p?1??j(p?1)(i?j)p]?O(hp?2)

(p?1)!j?0第11页

基于四次样条的矩阵微分方程数值解

于是局部误差的阶和方法的阶相同。

对于不同的p和i值，解方程(2.3.6)后，就得到不同标准的方法。 对p=2和i=1时，方程(2.3.6)变成 对于l=0，1+?1=0 对于l=1，1-?0-?1=0 对于l=2，1-2?0=0

由此得到?1=-1，?0=1/2，?1=1/2。这时对应的公式是

hyk?1=yk+[f(xk,yk)?f(xk?1,yk?1)] (2.3.9)

2这就是Euler法中的隐式格式，有时也称为修正的Euler公式。

重复应用公式(7.3.3),多步公式的局部误差和舍入误差都会有积累，可以证明，当且仅当方程 p(z)=?ajzi?j=0 (2.3.10)

j?0i在单位圆外无零点，在单位圆周上零点重数为1时，对于小的h值和小的舍入误差，则舍入误差和初值函数的误差函数的误差传播影响也是小的。具有这种性质的方法称为稳定的方法。 多步法的计算问题：

用k步法计算时，须知道初值y0,y1,...yk?1，其中y(x0)=y0是给定的初值，其余是附加初值。计算附加初值有几种方法。一种是单步法，比如Euler法和

Rung-Kutta法及其他单步法。为了保持多步法的精度，计算附加初值时要将

x0,xk?1之间节点加密或采用和多步法有同样阶的Rung-Kutta法。

多步法的第二个问题是如何选择阶p(或者步数k)。从收敛阶的观点，自然希望把p取大一些。但是高阶收敛方法要求解的光滑型也很高，否则达不到高精度的密度。并且高阶多步法的绝对稳定域也小，所以p的大小要考虑到解的光滑型和稳定性以及总的工作量。

多步法的第三个问题是步长h的选取。理论上似乎按照误差估计式选定h是合理的，但那种估计往往偏大，因此选定的h可能过小，既不必要也不经济。实

第12页

际用的步长h不是一次取定，而是根据精度要求，又粗到细逐渐调整（选步长），当h达到要求后就以此为步长计算。在计算中还可以改变步长，但计算过程变复杂了，这就是多步法的缺点。与此相反，单步法适合用变步长计算

2.4 外推法:

在求解初值问题

?y'?f(x,y) (2.4.1) ?y(x)?y00?的时候，可以写为

y(x)?y0??f(t,y?t?)dt (2.4.2)

0x而(2.4.2)式右端第二项的积分可以用各种各样的数值积分求解。这当然也可以用外推法来求（2.4.1）的解。

让我们回忆一下Richardson外推法。设T是某个被逼进的对象，T(h)是依赖于步长h的T的一个近似式 T?T(h) (2.4.3)

如果T与T(h)之间能表示为渐近式

T?T(h)+a1h?1?a2h?2?...?akh?k?... (2.4.4) 而0

T0(h)?T(h)?? (2.4.5) Tk?1(akh)?ak?kTk?1(h) k?1,2. ..??Tk(h)?1?ak?k?所以使用外推法的关键，是看存在什么样的渐近展开式(2.4.4)。对于Euler法，我们知道局部截断误差可写为

yk+1-y(yk+1 )=H( xk,h)= ??j?21(j)jyh (2.4.6) j!因此，经过几次计算得到T(h)后，鉴于全程截断函数误差比局部截断误差的阶

第13页

基于四次样条的矩阵微分方程数值解

数低一次，可设想y(x)?T(h)??ajhj (2.4.7)

j?1对于梯形法，因为它是由梯形求积公式推出来的，而梯形求积公式的截断误差具有展开式

?bhjj?12j?1

所以由前面的知识可知梯形法的局部截断误差为

y(yk?1)-xk?1??bjh2j?1 (2.4.8)

j?1如前，经过几次计算得到T(h)后，可设想

y(x)?T(h)??bjh2j (2.4.9)

j?1?同样，对于由中心商得到的Euler中心法公式 yk?1=yk?1+2h f(xk,yk) (2.4.10) (2.4.8)也成立。

根据(2.4.5)，令ak?(hi/hi?k)1/k，外推法计算公式为

T?Tiki?1k?1?1iTki?1?Tk?1 (2.4.11) ??(hi/hi?k)?1其中，对(2.4.7)，?k=k, ?=1,对(2.4.9)，?k=2k, ?=2,。注意T(h)是步长为h时求得的y(x)的近似值（例如yk?1为y(yk?1)在步长为h时的近似值，这时就取。与前面不同的是T0i=T(hi)。 T(h)=yk?1）

3 基于四次矩阵样条的微分方程求解算法

本文的工作是给出如(1)所示的一阶线性矩阵微分方程的一种数值求解方法

Y'(x) y'(x)= A(x)Y(x) + B(x), a < x < b,

第14页

Y(a)= Y a (1) 其中Ya,,Y?Cr?q,A:[a,b] ?Cr?r,B:[a,b] ?Cr?q满足A,B?C1([A,B]),它保证了（1）只有一个解，并且这个解是连续的、可微的 [10，P99]。

在参考文献[11]中，有许多不同的工程应用问题的数学模型都是以（1）的形式出现的。一些非线性的方程（如黎卡提方程[12-14]）通过线性化后也可以化为（1）的形式。利用固定步长的线性多步法来解决诸如（1）类型的方程的近似数值方法在参考文献[15]中有详细介绍，但是这个误差的范围依赖于步长h，且是h的若干次指数次幂，所以在实际操作中应使h尽可能小，并且在这些方法中为了得到一个连续的解，还必须要采用一些插值技巧。另外一些方法则是从Magnus or Fer方法发展起来的，利用这些方法求解问题（1）时，要计算一些矩阵的的幂，这就增加了计算成本。

文献[17]给出了用三次样条构造标量形式的常微分方程的数值解的方法，所得近似解具有以下优点：1、近似解在区间[a , b]上是C1连续的；2、近似解易于求值计算；3、近似解所达到的误差阶是O(h4)。最近，文献[18]也用了样条方法来构造另外一些标量型的微分方程的近似数值解。参考文献[19]研究了利用基于Hermite插值技巧的隐式样条方法来构造向量型的微分方程的数值解的问题。

在本文中，我们给出了利用四次矩阵样条来构造一阶矩阵微分方程（1）的近似数值解的方法。并与利用三次矩阵样条构造一阶矩阵线性微分方程的数值解的结果进行了比较，比较结果表明用四次矩阵样条所构造的数值解的逼近效果要比用三次矩阵样条所构造的数值解的逼近效果要好。所得方法仍然具有参考文献[17]中方法所具有的优点。

在这次工作中我们以Cp?q表示p?q阶复矩阵，如果A?Cr?s我们定义矩阵A的2次范数为A=supz?0Azz，在这里z是一个属于Cs的复向量,z=(ztz)1/2是z

的欧几里得范数，根据参考文献[20，p.56]，我们可得到A满足如下形式的不等式

maxijaij?A?srmaxijaij (2)

我们以Pn(x)表示以实变量x为自变量的n次矩阵多项式。设g是定义在区间[a,b]上

第15页

[0.7,0.8] ?1.00385?0.974971x?0.564182x2?0.0864059x3?0.0888587x4???0.0215663?0.860243x?1.35683x2?0.0573616x3?0.422203x4?????,?? ???4.12779?10?6 [0.8,0.9] [0.9,1] ??1.00709?0.95879x?0.594521x2?0.0611236x3?0.0967595x4?,??? ??2340.04076?0.764275x?1.53677x?0.0925888x?0.469063x??????5.10339?10?6 ??1.01471?0.9249x?0.651004x2?0.0192839x3?0.108382x4?,??? 6.23737?10?6 ??2340.0849726?0.567774x?1.86427x?0.335182x?0.53645x??????

表2

[xi,xi?1] 近似解 近似误差 [0,0.1] ???x23?1?x??0.172635x,????2?? ???x?x2?0.523192x3???????0.99977?1.00069x?0.493104x2?0.19562x3?,????? 23????0.0000967745?1.0029x?0.970968x?0.619967x???6.33769?10?6 [0.1,0.2] 6.34441?10?6 [0.2,0.3] ??0.999844?1.00269x?0.4831x2?0.212294x3?,?8.32925?10?6 ???? 23????0.000617244?1.01071x?0.931932x?0.685026x???[0.3,0.4] ??0.99916?1.00953x?0.460305x2?0.237623x3?,????? 23????0.00376114?1.04215x?0.827136x?0.801466x???8.34117?10?6 第21页

基于四次样条的矩阵微分方程数值解

[0.4,0.5] ??0.997697?1.0205x?0.432875x2?0.260481x3?,?11.5366?10?6 ???? 23????0.00954449?1.08552x?0.718698x?0.891831x???[0.5,0.6] ??0.994112?1.04201x?0.389854x2?0.289161x3?,?11.5579?10?6 ???? 23????0.0273079?1.19211x?0.505537x?1.03394x???[0.6,0.7] ??0.987639?1.07438x?0.335913x2?0.319129x3?,????? 23????0.0536042?1.32359x?0.286401x?1.15568x???16.368?10?6 [0.7,0.8] ??0.976268?1.12311x?0.266294x2?0.35228x3?,????? 23????0.113731?1.58127x?0.0817239x?1.33098x???16.4047?10?6 [0.8,0.9] ??0.956625?1.19677x?0.174217x2?0.390646x3?,????? 23????0.196049?1.88996x?0.467586x?1.49175x???23.3063?10?6 [0.9,1] ??0.928364?1.29097x?0.0695481x2?0.429412x3?,?24.7135?10?6 ???? 23?????0.354975?2041972x?1.0562x?1.70976x??比较三次样条和四次样条的计算误差，可以看出四次样条的计算误差是三次样条的十分之一，计算精度比三次样条提高许多。对于相同的数值逼近问题，如果取相同的步长，四次样条的精度比三次样条高一个数量级，也就意味着如果要取得相同的计算精度，四次样条的计算次数只是三次的一半，所以利用四次样条可以减少计算量。 例4.2对于矩阵问题

?A(x)Y(?x)B ( Y'(x) x

第22页

?11???Y(0)??0?1?,x?[0,1]

?00?????1?x0?1?ex?x????x1? 其中 A(x)??exx?0??1e????2?(?3?ex)x?1?x?x2?ex(2?x)???B(x)??x?ex(1?x)e2x?xex?(5?x)(1?x2)?

??1?xex?1?x(5?x)??

并且该问题的精确解是

?1?xex?x???Y(x)??0?1?5x?x2?

?x?0??因为maxx?[0,1]||A2(x)?A'(x)||?5所以我们取M=5，由此得到下式

Y''(x)?[A'(x)?A2(x)]Y(x)?A(x)B(x)?B'(x)

Y'''(x)=[ A(x)A'(x)+A'(x)A(x)+A''(x)] Y(x)+[A'(x)+A2(x)]Y'(x)+ A(x)B'(x)

+ A'(x)B(x)+ B''(x)

?12???Y'(0)?A(0)Y(0)?B(0)??05?，

?10???并且可得

?01??01?????Y''(0)??02? ，Y'''(0)??00?

?00??00?????对于n=10,它满足n?M(b?a)/4，所以h?(b?a)/n?0.1。

利用前面的算法同样可以得到用三次样条和四次样条逼近的对比结果，其中四次样条逼近的结果如表3所示，三次样条的逼近结果如表4所示。

表3

第23页

基于四次样条的矩阵微分方程数值解

[xi,xi?1] 近似解 近似误差 [0，0.1] ???x2x3?1544?1?x?3.70075?10x,1?2x???0.430911x,????26???????1542?74???3.69393?10x,?1?5x?x?1.06649?10x?,? ???184?74x?6.50969?10x,?5.31896?10x????????[0.1,0.2]??1?x?8.88059?10?16x2?7.505.92039?10?15x3?1.11002?10?14x4,????5.76994?10?8 ????,??1?1.99998x?0.50036x2?0.16427x3?0.490835x4?????????13?17?162?153?154????1.10764?10?4.43058?10x?6.64987?10x?4.43098?10x?7.38251?10x,?? ???,?2?83?74?????1?5x?x?2.77673?10x?1.76067?10x?????162?163?154?????0.?1.x?1.1263?10x?7.50854?10x?1.87062?10x,???????8.95311?10?11?3058124?10?9x?5.37187?10?8x2?3.58124?10?7x3?3.63415?10?7x4??????1.20151?10?7 [0.2,0.3] ??1?1.x?4.43668?10?15x2?1.18287?10?14x3?1.10862?10?14x4,??????,????1.00001?1.99988x?0.501248x2?0.161308x3?0.0527852x4? ????????17?4.2769?10?16x?2.87538?10?15x2?7.36932?10?15x3?7.36736?10?15x4,?????2.24922?10?????,?2?83?74?????1?5.x?1.x?7.97667?10x?2.41066?10x?????17?1.x?2.56534?10?15x2?8.17571?10?15x3?9.28759?10?15x4,????2.77556?10?????????1.44123?10?9?2.7034?10?8x?1.75896?10?7x2?4.07257?10?7x3?5.93313?10?7x4??????1.88255?10?7 [0.3,0.4] ??1?1.x?1.51916?10?15x2?1.40646?10?15x3?5.6847?10?17x4,??????,????1.00006?1.9991x?0.505034x2?0.152894x3?0.0597971x4? ????????16?2.25912?10?15x?1.05587?10?14x2?2.24841?10?14x3?1.75105?10?14x4,????1.79019?10?????,?2?73?74?????1?5.x?1.x?2.38049?10x?3.72968?10x?????16?142?143?144?????1.x?1.34747?10x?2.74688?10x?2.04162?10x,??3.33067?10???????5.53712?10?9?6.60107?10?8x?2.89328?10?7x2?6.26572?10?7x3?2.68212?10?7x4??????2.63342?10?7 第24页

?142?143?144??1?x?1.23648?10x?1.94825?10x?1.12407?10x,??[0.4,0.5]?,? ?3.47323?10?7 ??234??1.00019?1.99786x?0.50967x?0.145168x?0.0646256x?????????9.76387?10?16?9.29493?10?15x?3.2769?10?14x2?4.97287?10?14x3?2.76225?10?14x4,????? ????2?73?74??1?5x?x?5.05176?10x?5.39923?10x????????15?142?143?144??1.38778?10?x?4.71066?10x?7.35?10x?4.26893?10x,?????8?7?72?73?74?????2.03462?10?1.92822?10x?6.81296?10x?9.91133?10x?7.42854?10x??????[0.5,0.6] ??1?x?4.7952?10?15x2?3.39747?10?15x3?1.99306?10?16x4,????????,??1.0007?1.99374x?0.522056x2?0.128654x3?0.0728829x4?????? ???15?14?142?143?144????3.04806?10?2.29006?10x?6.38176?10x?7.90535?10x?3.676896?10x,??????2?1.1034?10?6x3?8.39034?10?7x4????1?5x?0.999999x??????15?132?133?144??????4.996?10?x?1.06109?10x?1.30787?10x?5.94543?10x,???????2.79883?10?8?1.93854?10?7x?4.78732?10?7x2?5.5557?10?7x3?3.04978?10?8x4??????4.23943?10?7 [0.6,0.7] ??1?x?3.58159?10?14x2?3.78649?10?14x3?1.45607?10?14x4,??????,????1.00151?1.98537x?0.535458x2?0.113763x3?0.0790876x4? ????????15?14?142?133?144???4.21255?10x?9.87477?10x?1.01575?10x?3.84931?10x,????6.70586?10?????2?63?64?????1?5x?0.999999x?2.16636?10x?1.28193?10x?????14?x?1.75732?10?13x2?1.82369?10?13x3?7.10277?10?14x4,????1.16573?10?????????1.21871?10?7?8.05209?10?7x?2.01893?10?6x2?2.2196?10?6x3?1.12593?10?6x4??????5.5423?10?7 [0.7,0.8] ??1?x?5.94321?10?15x2?1.90568?10?15x3?3.56965?10?16x4,????????,??1.00385?1.97496x?0.564206x2?0.0863837x3?0.0888658x4?? ???????14?14?132?133?144???7.09874?10x?1.43637?10x?1.29268?10x?4.39507?10x,???1.30889?10?????2?63?64?????1?5x?0.999996x?4.29102?10x?2.04074?10x?????14?132?133?144?????x?2.76579?10x?48403?10x?8.28195?10x,??2.5091?10???????5.20669?10?9?7.90503?10?8x?4.62871?10?7x2?7.37648?10?7x3?5.96555?10?7x4??????6.873?10?7 [0.8,0.9] ??1?x?6.60126?10?14x2?5.19635?10?14x3?1.52861?10?14x4,??????,????1.00708?1.95881x?0.594487x2?0.0611496x3?0.0967514x4?? ???????14?13?132?133?144???1.01258?10x?1.79322?10x?1.39865?10x?4.01534?10x,????2.13601?10?????2?63?64?????1?5x?0.999991x?8.19135?10x?3.2596?10x?????14?x?3.99038?10?13x2?3.14611?10?13x3?9.31222?10?14x4,????4.70735?10?????????6.40389?10?7?3.99038?10?6x?6.51533?10?6x2?5.78137?10?6x3?2.17272?10?7x4??????8.51812?10?7 第25页

基于四次样条的矩阵微分方程数值解

[0.9,1] ??1?x?8.29596?10?14x2?5.83863?10?14x3?1.53666?10?14x4,????????,??1.01472?1.92488x?0.651037x2?0.0192607x3?0.108387x4?????? ???14?5.68534?10?14x?8.4196?10?14x2?5.53339?10?14x3?1.40685?10?14x4,??1.42149?10????????23?64?????1?5.00001x?0.999981x?0.000015782x?5.35038?10x?????14?132?133?134?????8.2606?10?x?5.60913?10x?3.96464?10x?1.04398?10x,?????????7.54631?10?7?3.92561?10?6x?7.54631?10?6x2?6.54505?10?6x3?2.38485?10?7x4??????1.0636?10?6 表4

[xi,xi?1] [0，0.1] ??1?x?0.x3,近似解 近似误差 ??????2,???1?2x?x?0.172313x3??????2???3??0.x,???? ,????23??1?5x?x?0.0000511502x???????3??x?0.x,???????63???1.77022?10x????????1.39627?10?6 [0.1,0.2] ?1?1.x,???????,??0.999976?2.00073x?0.492747x2?0.19649x3???????? ???18?2.03067?10?16x?2.03067?10?15x2?6.76889?10?15x3,?????6.76889?10?????,?23?????1?4.9999x?1.00006x?0.000158467x?????172?163????x?2.35206?10x,?0.?1.x?7.05617?10??????10?9?82?63??1.54825?10??4.64475?10x?4.64475?10x?1.92505?10x?????1.39627?10?6 [0.2,0.3] ??1?1.x?2.13026?10?14x2?2.8386?10?14x3,????????,??0.999861?2.00245x?0.484122x2?0.210864x3?????? ???16?2.21437?10?15x?1.00565?10?14x2?1.33764?10?14x3,????1.54394?10?????,?23?????1?5.00005x?0.999799x?0.000281181x?????162?163????0.?1.x?3.50695?10?x?4.66888?10x,????????7.98671?10?8?1.19568?10?6x?5.9552?10?6x2?8.0777?10?6x3??????1.43684?10?6 第26页

[0.3,0.4] ??1?1.x?2.98265?10?14x2?2.84241?10?14x3,????????,??0.999085?2.01021x?0.458261x2?0.239599x3?????? ???16?8.86587?10?15x?2.68776?10?14x2?2.76615?10?14x3,?????9.5363?10?????,?23??????0.999984?4.99986x?1.00044x?0.000425581x????16?152?153?????1.x?4.44774?10x?4.8647?10x,?1.11022?10??????7?623??6.6746?10?6.27759?10x?0.0000189557x?0.0000196011x??????1.43614?10?6 [0.4,0.5] ??1?1.x?3.82909?10?14x2?2.83403?10?14x3,????????,??0.997911?2.01901x?0.436255x2?0.257937x3?????? ???15?3.15699?10?14x?7.42117?10?14x2?5.65796?10?14x3,??4.4378?10???????,?23?????1.00005?5.00034x?0.999213x?0.000593048x?????16?152?153?????1.x?6.13427?10x?3.95364?10x,??4.996?10????????2.73762?10?6?0.0000192605x?0.0000448896x2?0.0000336034x3????????1?1.x?3.75051?10?14x2?2.78165?10?14x3,????????,??0.99362?2.04476x?0.384763x2?0.292265x3?????? ???15?1.19489?10?14x?1.28259?10?14x2?1.44546?10?15x3,?????2.81533?10?????,?23?????0.999876?4.99931x?1.00129x?0.000792557x?????16?152?153?????1.x?9.47701?10x?6.45388?10x,?7.21645?10??????623??8.35677?10?0.0000473058x?0.0000882431x?0.0000551518x??????1.49859?10?6 [0.5,0.6]1.49774?10?6 [0.6,0.7] ??1?1.x?2.63842?10?13x2?1.39598?10?13x3,????????,??0.98862?2.06976x?.034309x2?0.315417x3????????? ?14?14?132?143??9.45552?10x?1.50503?10x?7.79327?10x,?????1.93366?10???,?23?????1.00027?5.00127x?0.998018x?0.00102612x?????15?142?153?????2.77556?10??1.x?1.85495?10x?9.11641?10x,????????0.000020469?0.0000968231x?0.000151972x2?0.0000783009x3????????1?1.x?4.39029?10?13x2?1.95102?10?13x3,????????,??0.97449?2.13032x?0.256581x2?0.356611x3?????? ???14?8.43694?10?14x?1.05104?10?13x2?4.37846?10?14x3,????2.24125?10?????,?23?????0.999468?4.99784x?1.00292x?0.00130831x?????15?142?153?????1.x?1.4073?10x?6.41813?10x,?2.55351?10???????0.0000451908?0.000184576x?0.000250027x2?0.000113127x3??????1.57321?10?6 [0.7,0.8]1.57205?10?6 第27页

基于四次样条的矩阵微分方程数值解

[0.8,0.9] ??1?1.x?4.28079?10?13x2?1.66193?10?13x3,????????,??0.959648?2.18597x?0.187012x2?0.385599x3?????? ???14?1.90541?10?13x?2.37818?10?13x2?9.90823?10?14x3,????5.07249?10?????,?23?????1.00098?5.0035x?0.995844x?0.00163979x?????15?142?143?????1.x?2.9315?10x?1.16602?10x,??6.77236?10????????0.0000876052?0.000313409x?0.000372454x2?0.00014624x3????????1?1.x?3.89588?10?13x2?1.36647?10?13x3,????????,??0.923466?2.30658x?0.0530041x2?0.435231x3?????? ???13?6.34353?10?13x?6.78731?10?13x2?2.4038?10?13x3,?????1.96743?10?????,?23?????0.998289?4.99454x?1.0058x?0.00204807x?????14?142?143?????1.x?6.99077?10x?2.50889?10x,?2.0095?10???????0.000164691?0.000527578x?0.000561976x2?0.000199845x3??????1.66154?10?6 [0.9,1] 1.665985?10?6 比较上面两表中数据，同样可以看出用四次样条比用三次样条的计算精度高，四次样条的逼近效果更好。

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 Nostra~d, Reinhold, (1971). 3. W.T. Reid, Riccat~ Differential Equations, Academic Press, (1972). 4. R.A. LaBudde, Classical Mechanics of Molecular Collisions, University of

第28页

Wisconsin Theretica] Chemistry Inst. Report WIS-TC1-414, (1973).

5. M. Scott, In~ariant ~mbeddiz~g and i~s Applications to Ordinary Differential Equations, Addison-Wesley,(1973).

6. P. Marzulli, Global error estimates for the standard parallel shooting method, J. Computat. Appl. Ma~hs.34, 233-241, (1991).

7. K. Raktorys, The Method of Di~cre$~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 Produc~s 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 E~ations, 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. Cort~s, Rational matrix approximation with a priori error bounds for non-symmetric matrix Riccati Equations with analytic coefficients, IMA J. Numer. Anal. 18 (4), 545-561, (1998).

15. L. Jddar and E. Ponsoda, Continuous numerical solutions and error bounds for matrix differential equations,[n Int. Prec. b-b/st Int. Coil. Num. Anal., (Edited by D. Bainov and V. Covachev), pp. 73--88, VSP. Utrech,The Netherlandsj (1993). 16. S. Blanes, F. Casas, J.A. Oteo and J. Ros, Magnus and Fer expansion for matrix differential equations: the convergence problem, J. Phys. AppL 31,259--268, (19.q8). 17. F.I~ Loscalzo and T.D. Talbot, Spline function approximations for solutions of

第29页

基于四次样条的矩阵微分方程数值解

ordinary differential equations, SIAM J. Numer. Anal. 4 (3), 433-445, (1967). 18. E.A. Al-Said and M.A. Noor, Cubic splines method for a system of third-order boundary value problems, AppL Math. Computation 142, 195--204, (2003). 19. Gh. Micula and A. Revnlo, An implicit numerical spline method for systems for ODE's, Appl. Math. ComputationII15 121-132, (2000).

20- G.H. Golub and C.F. VanLoan, Ma~riz Computations, Third Edition, Johns Hopkins University Press,(1996).

21. E. Dcfez, A. Hervtts, A. Law, J. Villanueva-OIler and R. Villanueva, Matrix cubic splines for progressive transmission of images, Journal of Ma~hemahcal Imaging and VLs~on 17 (1), 41-53, (2002).

第30页

本文来源：https://www.bwwdw.com/article/lvu5.html