常微分方程初值解问题数值解的实现和分析

1、 课程设计说明书（论文） 第 PAGE 46页 常微分方程初值解问题数值解的实现和分析摘要此次课程设计是简单介绍几种常微分方程的初值问题数值解的求法以及其matlab实现。在这次课程设计中我对几种常见的数值解法做了简单介绍如Eulor方法，改进Eulor方法，RungeKutta方法，Adams方法等，在这里我们仅给予了简单的推到和介绍。对于一些方法我给出相应的流程图以及简单的算法分析，最后参考多种方法选取了三种典型的代表方法并给出相应的程序。最后我们根据具体的要求针对方程： 我给出相应的程序编写和使用方法，仅供参考使用。关键词：微分方程的数值解，Eulor方法，RungeKutta方法，Ad

2、ams预测校验目 录 TOC o 1-3 h z u HYPERLINK l _Toc295630417 1 引言 PAGEREF _Toc295630417 h 1 HYPERLINK l _Toc295630418 2 常用的方法 PAGEREF _Toc295630418 h 2 HYPERLINK l _Toc295630419 2.1 Euler方法 PAGEREF _Toc295630419 h 2 HYPERLINK l _Toc295630420 2.2 改进的Euler方法 PAGEREF _Toc295630420 h 5 HYPERLINK l _Toc295630421

3、 2.3 Runge-Kutta方法 PAGEREF _Toc295630421 h 5 HYPERLINK l _Toc295630422 2.3.1 二级Runge-Kutta法 PAGEREF _Toc295630422 h 6 HYPERLINK l _Toc295630423 2.3.2 常用四阶Runge-Kutta 法 PAGEREF _Toc295630423 h 6 HYPERLINK l _Toc295630424 2.4 线性多步法 PAGEREF _Toc295630424 h 7 HYPERLINK l _Toc295630425 2.4.1 Adams方法 PAGE

4、REF _Toc295630425 h 8 HYPERLINK l _Toc295630426 2.4.2 Adams预测校正系统 PAGEREF _Toc295630426 h 10 HYPERLINK l _Toc295630427 2.5 一阶微分方程组和高阶微分方程 PAGEREF _Toc295630427 h 10 HYPERLINK l _Toc295630428 3 程序开发思路及简单注解 PAGEREF _Toc295630428 h 11 HYPERLINK l _Toc295630429 3.1 改进欧拉算法级程序 PAGEREF _Toc295630429 h 12 H

5、YPERLINK l _Toc295630430 3.2 四阶龙格库塔算法及其程序： PAGEREF _Toc295630430 h 13 HYPERLINK l _Toc295630431 3.3 Adams预测校正系统 PAGEREF _Toc295630431 h 16 HYPERLINK l _Toc295630432 4 求解例题 PAGEREF _Toc295630432 h 19 HYPERLINK l _Toc295630433 4.1 高阶方程的转化和真实解的求解 PAGEREF _Toc295630433 h 19 HYPERLINK l _Toc295630434 4.2

6、 改进Eulor方法求解 PAGEREF _Toc295630434 h 20 HYPERLINK l _Toc295630435 4.3 四阶RungeKutta方法求解 PAGEREF _Toc295630435 h 21 HYPERLINK l _Toc295630436 4.4 四阶Adams预测校正系统方法求解 PAGEREF _Toc295630436 h 23 HYPERLINK l _Toc295630437 5 结果及误差分析 PAGEREF _Toc295630437 h 24 HYPERLINK l _Toc295630438 小结 PAGEREF _Toc2956304

7、38 h 26 HYPERLINK l _Toc295630439 参考文献 PAGEREF _Toc295630439 h 27 HYPERLINK l _Toc295630440 附录： PAGEREF _Toc295630440 h 28 HYPERLINK l _Toc295630441 附录1 PAGEREF _Toc295630441 h 28 HYPERLINK l _Toc295630442 附录2 PAGEREF _Toc295630442 h 29 HYPERLINK l _Toc295630443 附录3 PAGEREF _Toc295630443 h 30 HYPERL

8、INK l _Toc295630444 附录4 PAGEREF _Toc295630444 h 31 HYPERLINK l _Toc295630445 翻译 PAGEREF _Toc295630445 h 33 1 引言微分方程数值解一般可分为：常微分方程数值解和偏微分方程数值解。自然界与工程技术中的许多现象，其数学表达式可归结为常微分方程（组）的定解问题。一些偏微分方程问题也可以转化为常微分方程问题来（近似）求解。Newton最早采用数学方法研究二体问题，其中需要求解的运动方程就是常微分方程。许多著名的数学家，如 Bernoulli（家族），Euler、Gauss、Lagrange和Lap

9、lace等，都遵循历史传统，研究重要的力学问题的数学模型，在这些问题中，许多是常微分方程的求解。为什么要研究数值解法？我们知道，只有少数十分简单的微分方程能够用初等方法求得它们的解，多数情形只能利用近似方法求解。此外自然科学各学科领域中的大量现象和问题，其其数学的表达式归结为微分方程的定解问题。如何通过如何数值方法获得微分方程问题的满足指定精度的近似解，是科学与工程中的一类最基本、应用最广泛的计算问题。微分方程的初值问题是微分方程定解问题的一个典型代表。常微分方程初值问题中最简单的例子是人口模型，设某一地区在时的人口为为已知的，该地区的人口自然增长率为，人口的增长与人口总数成正比，所以t时刻的

10、人口总数为，满足以下微分方程：很多物理系统与时间有关，从卫星运行轨道到单摆运动。从化学反应到物种竞争都是随着 时间的连续而不断变化的。微分方程是描述连续变化的数学语言，微分方程的就是确定满足给定方程的可微函数，研究他的数值方法。在常微分方程课中已经讲过的级数解法，逐步逼近法等就是近似解法。这些方法可以给出解的近似表达式，通常称为近似解析方法。还有一类近似方法称为数值方法，它可以给出解在一些离散点上的近似值。利用计算机解微分方程主要使用数值方法。我们考虑一阶常微分方程初值问题：在区间a, b上的解，其中f (x, y)为x, y的已知函数，y0为给定的初始值，将上述问题的精确解记为y(x)。数值

11、方法的基本思想是：在解的存在区间上取n + 1个节点这里差，i = 0,1, , n称为由xi到xi+1的步长。这些hi可以不相等，但一般取成相等的，这时。在这些节点上采用离散化方法，（通常用数值积分、微分。泰勒展开等）将上述初值问题化成关于离散变量的相应问题。把这个相应问题的解yn作为y(xn)的近似值。这样求得的yn就是上述初值问题在节点xn上的数值解。一般说来，不同的离散化导致不同的方法。2 常用的方法在学习常微分方程数值解时我们会遇见不同情况的解法，在这里我们仅就几种简单的、常见的数值积分法做简单是介绍。2.1 Euler方法考虑初值问题为了求得它在等距离散点上的数值解，首先将(2.1

12、)离散化。设，将式(2.1)离散化的办法有Taylor展开法、数值微分法及数值积分法。如果在点将作Taylor展开,得 (2.3)那么当充分小时，略去误差项，用近似替代、近似替代，并注意到，便得： (2.4)其中. 用(2.4)求解(2.1)的方法称为Euler方法。如果利用差商近似替代微商，那么可得： (2.5)在(2.5)中若用近似替代、近似替代，同样得到递推公式(2.4)。如果在上对积分，得： (2.6)那么对上式右端积分用左矩形求积公式，并用近似替代、近似替代，也可得到递推公式(2.4)。我们知道，在平面上，微分方程的解称为积分曲线，积分曲线上一点的切线斜率等于函数的值。如果在D中每一

13、点，都画上一条以在这点的值为斜率并指向增加方向的有向线段，即在D上作出了一个由方程确定的方向场，那么方程的解. 从几何上看,就是位于此方向场中的曲线，它在所经过的每一点都与方向场的该点的方向相一致。从初始点出发，过这点的积分曲线为，斜率为. 设在附近可用过点的切线近似表示，切线方程为. 当时，的近似值为，并记为，这就是得到时计算的近似公式当时，的近似值为，并记为. 于是就得到当时计算的近似公式重复上面方法，一般可得当的计算的近似公式：如果，则上面公式就是(2.4). 将连续起来，就得到一条折线，所以Euler方法又称为折线法。由公式(2.4)看出，已知便可算出. 已知，便可算出，如此继续下去，

14、这种只用前一步的值便可计算出的递推公式称为单步法。若在(2.6)中，其右边的积分由数值积分的右矩形公式近似，并用替代，替代，则可得到： (2.7)并称(2.7)为后退Euler公式。Euler公式(2.4)是关于的一个直接计算公式，是显式的；而公式(2.7)右端是含有的一个函数方程，因此是隐式的，也是单步法。图1若在公式(2.6)中,其右边积分用数值梯形积分公式近似，并用替代，替代，则可得到梯形方法公式 (2.8)梯形方法同后退Euler方法一样都是隐式单步法.对于隐式方法，通常采用迭代法1。对后退Euler方法，先Euler方法计算，并将它作为初值，即,再把它代入(1.7)的右端，便得到后退

15、Euler方法的迭代公式为 (2.9)同样地，仍用Euler方法提供初始值，梯形方法的迭代公式为 (2.10)2.2 改进的Euler方法梯形方法的迭代公式(2.10)比Euler方法精度高,但其计算较复杂,在应用公式(2.10)进行计算时,每迭代一次,都要重新计算函数的值,且还要判断何时可以终止或转下一步计算.为了控制计算量和简化计算法,通常只迭代一次就转入下一步计算.具体地说,我们先用Euler公式求得一个初步的近似值,称之为预测值,然后用公式(2.10)作一次迭代得,即将校正一次.这样建立的预测校正方法称为改进的Euler方法:预测: 校正: （2.10）2.3 Runge-Kutta方

16、法我们利用在点的Taylor展开的前两项导出Euler公式. 很自然地，我们想到：若将在点的Taylor展开多取几项，则有希望获得高阶方法。但是直接利用Taylor展开取多项的办法，需要计算的高阶导数，计算量较大。因此在建立Runge-Kutta方法时不采用求高阶导数的办法，而是用计算不同点上的函数值，然后对这些值作线性组合,构造近似公式，把近似公式和解的Taylor展开相比较，使它们在前面的若干项相同，从而使近似公式达到一定的阶2。采用区间内若干点对这些值作线性组合，一般形式为 (2.11)其中为待定参数，依次类推方法称为显型r级Runge-Kutta方法。2.3.1 二级Runge-Kut

17、ta法当r=2时，将f(x,y)在(xk , yk)上展开 (2.12)又因 (2.13) 以上两式代入2.13又由得：这里三个方程四个未知量。因此，具有无穷多个数量组。特殊的若取 相应的二级公式为：显然公式等价于： (2.14)此公式称为中点法。2.3.2 常用四阶Runge-Kutta 法通常人们所说的龙格库塔法是指四阶而言的。我们可以仿照二阶的情形推导出此公式，不过太繁杂，此处从略，常用的四阶公式是 (2.15) 此外常用的四阶Runge-Kutta公式还有Gill方法：(2.16)2.4 线性多步法前面介绍的方法，统称为单步法，就是在计算yn+1时，只用到前面一步yn的值。这里我们将介

18、绍利用前边已经算出来的若干个值yn-k，yn-k+1，yn-1，yn，来求得yn+1的高精度公式线性多步方法。我们已经知道初值问题 (2.1)与积分方程 (2.17)等价。下面的方法，其基本思想是用一个插值多项式P(x)来代替（2.17）中的被积函数f (x, y)，然后用 (2.18)代替(2.17)。这种方法实际上要分两步来做：1求出开头几个点上的近似值，即计算“表头”；2利用(2.18)逐步求后面点上的值。选取不同的插值点作插值多项式，就会得出不同的数值解法。2.4.1 Adams方法我们这里只讨论其中的一种，叫阿当姆斯（Adams）方法。 Adams方法用等步长来数值求解（1.1），(

19、1.2)，为正数。设已用某种方法求出来了。已是知道。步显示Adams方法：用经过点的插值多项式(2.19)来近似并取我们就建立了步显式Adams公式 (2.20)其中 (2.21) 我们将几个低阶显式Adams方法的公式及主局部截断误差的系数例在下表2中。步数方法阶公 式12341234表1步隐式Adams方法：用经过点的插值多项式 (2.22) 来近似在进行数值积分，可以得到步隐式Adams方法 (2.23)其中系数由下述积分给出 (2.24) 我们在此给出几个个低阶隐式Adams方法的公式及主局部截断误差的系数列在表2-3中。步数方法阶公 式12233445表22.4.2 Adams预测校

20、正系统上述给出的Adams显式方法计算简单，但精度比隐式方法差，而隐式方法由于每步要做迭代，计算不方便.为了避免迭代，通常可将同阶的显式Adams方法与隐式Adams方法结合，组成预测-校正方法.以四阶方法为例，可用显式方法(7.5.8)计算初始近似，这个步骤称为预测(Predictor)，以P表示，接着计算f值(Evaluation),，这个步骤用E表示，然后用隐式公式(2.23)计算，称为校正(Corrector)，以C表示，最后再计算，为下一步计算做准备.整个算法如下： 预测 P：求值 E：校正 C：求值 E： (2.24)公式称(2.24)为四阶Adams预测-校正方法(PECE)3。

21、2.5 一阶微分方程组和高阶微分方程所有应用于一阶常微分方程初值问题的数值方法都可以直接推广到一阶常微分方程组，只需把公式中的未知函数改为向量形式。高阶方程总可以降为一阶方程组，进而 用数值方法求解。一般的一阶常微分方程组初值问题是求解 （2.25）（2.25）的向量形式是：其中。高阶常微分方程初值问题一般为： (2.26)其中是给定多元函数，为给定值。引进新的变量函数： (2.27)则初值问题（2.26）化成了一阶常微分方程组初值问题 (2.28) 通过求解(2.28) 得到（2.26）的解。3 程序开发思路及简单注解介于本次课程设计求解常微分方程数值解所涉及的方法比较多，这里我们进选取有代

22、表的三中方法做出程序开放思路与注解，其它均可参照自行开发。3.1 改进欧拉算法级程序算法公式： 算法流程图：开始输入h和区间确定迭代次数n计算计算迭代次数kkn输出k=n图2算法分析：由于我们进行的是固定步长h的迭代运算所以我们要先计算出迭代次数；其次进入循环主体：每次循环我们需要求出，然后再求出，最后求出；最后我们当达到迭代次数时我们退出循环体（方程组类似此处不赘余）。在此我们给出方程组的求解部分主体循环程序：for k=1:n %进入循环体f1=h*feval(dydx,X(k),Y(k,:); %求出f1=f1;x2=X(k);y2=Y(k,:);f2=h*feval(dydx,x2,Y

23、(k,:)+f1); %求出f2=f2;Y(k+1,:)=Y(k,:)+(f1+f2)/2 % end退出循环体。3.2 四阶龙格库塔算法及其程序：四阶Runge-Kutta的计算公式：对于方程组我们可以将其向量化，例如对于二阶方程：引进新的变量，即可化为下列一阶方程组的初值问题：针对这个问题应用四阶龙格库塔公式，有对应的分量形式为：算法分析及其流程图(向量图类似)；开始输入h和区间确定迭代次数n计算计算迭代次数kkn输出k=n计算计算图3算法分析：由图3.2程序开放时我们的运算过程如下： 计算需要循环的次数；依此计算： 计算：如果达到循环次数没有达到则返回（2）如果达到则停止计算。对于方程组

24、同样是如此循环计算只是每一次计算值多余现在的在此我们不再重复叙述。微分方程组的数值解部分程序开发及其分析：写成分量形式：for i=2:n K(1,1)=f1(x(i-1),y(i-1),z(i-1); %计算 K(2,1)=f2(x(i-1),y(i-1),z(i-1); %计算 K(1,2)=f1(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); %计算 K(2,2)=f2(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); %计算 K(1,3)=f1(x(i-1)+h/2,y(i-1)+K(1,

25、2)*h/2,z(i-1)+K(2,2)*h/2); %计算 K(2,3)=f2(x(i-1)+h/2,y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2); %计算 K(1,4)=f1(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); %计算 K(2,4)=f2(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); %计算 y(i)=y(i-1)+h/6*(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4); %计算 z(i)=z(i-1)+h/6*(K(2,1)+2*K(2,2)+2*K(2

26、,3)+K(2,4); %计算end编辑成为向量形式：for k=1:nk1=feval(dydx,X(k),Y(k,:); x2=X(k)+h/2;y2=Y(k,:)+k1*h/2;k2=feval(dydx,x2,y2); k3=feval(dydx,x2,Y(k,:)+k2*h/2); k4=feval(dydx, X(k)+h,Y(k,:)+k3*h); Y(k+1,:)=Y(k,:)+h*(k1+2*k2+2*k3+k4)/6;k=k+1;end运行时每一句与分量形式相同在此不做解释了。3.3 Adams预测校正系统Adams预测校正系统的计算公式：预测 P：求值 E：校正 C：求值

27、 E：对与方程组我们只需要像四阶龙格库塔公式那样将其改为向量或分量形势即可这里不做重复。 算法分析：由下图的流程图我们可以得出在进行程序的编写需要进行如下步骤：进行表头的计算；求出预测值P；求出值E；求出校正值C；在此求值E；如果达到循环次数则结束，如果达不到则返回(2)。程序流程图： 开始输入h和区间确定迭代次数n计算预测值kn输出k=n计算值E计算计算表头计算校正值求值迭代次数k图4Adams预测校正系统的部分程序(向量形式)4：for k=1:3k1=feval(dydx,X(k),Y(k,:); x2=X(k)+h/2;y2=Y(k,:)+k1*h/2;k2=feval(dydx,x2

28、,y2); k3=feval(dydx,x2,Y(k,:)+k2*h/2); k4=feval(dydx,X(k)+h,Y(k,:)+k3*h); Y(k+1,:)=Y(k,:)+h*(k1+2*k2+2*k3+k4)/6;end%采用四阶龙格库塔法计算表头。X;Y;for i=1:4f(i,:)=feval(dydx,X(i),Y(i,:);endfor k=4:nP=Y(k,:)+(h/24)*(-9)*f(k-3,:)+37*f(k-2,:)-59*f(k-1,:)+55*f(k,:);%计算预测值f(k+1,:)=feval(dydx,X(i+1),P);%计算函数值Y(k+1,:)=

29、Y(k,:)+(h/24)*(f(k-2,:)-5*f(k-1,:)+19*f(k,:)+9*f(k+1,:); %求校正值f(k+1,:)= feval(dydx,X(k+1),Y(k+1,:);%求出最后值end4 求解例题现有某一周期外力作用的共振弹簧系统的微分方程为： 4.1将他转化为一阶方程组后采用：改进Eulor方法，四阶RungeKutta方法，和三步四阶Adams预测校正系统分别求解在x=0：0.05：1上关于位移y和速度的数值解并画出位移和速度的图像。4.1 高阶方程的转化和真实解的求解我们将其转化为一阶方程组的形势为了便于观察函数的真实解我们可以首先求出方程的真实解然后与数

30、值解对比。真实解的求解：在matlab输入窗口输入：y=dsolve(D2y+64*y-16*cos(8*x)=0,y(0)=0,Dy(0)=0,x)输出：y =sin(8*x)*x4.2 改进Eulor方法求解编写改进Eulor函数文件：gjeuler保存；建立函数文件dydx并保存；编写M文件gaijineuler并保存。调用gaijineuler 文件运行并输出结果：k = 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21X =0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0

31、.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9500 1.0000Y = 0 0 0.0200 0.8000 0.0768 1.4088 0.1551 1.6077 0.2303 1.2728 0.2749 0.4107 0.2651 -0.8348 0.1874 -2.2063 0.0433 -3.3834 -0.1493 -4.0498 -0.3578 -3.9654 -0.5405 -3.0262 -0.6547 -1.3003 -0.6656 0.9689 -0.5546 3

32、.3962 -0.3249 5.5195 -0.0037 6.8857 0.3607 7.1412 0.7063 6.1110 0.9675 3.8486 1.0876 0.6456P = 1.0000 0 0 0 0 2.0000 0.0500 0.0200 0.8000 0.0200 3.0000 0.1000 0.0768 1.4088 0.0568 4.0000 0.1500 0.1551 1.6077 0.0782 5.0000 0.2000 0.2303 1.2728 0.0752 6.0000 0.2500 0.2749 0.4107 0.0446 7.0000 0.3000 0

33、.2651 -0.8348 0.0098 8.0000 0.3500 0.1874 -2.2063 0.0777 9.0000 0.4000 0.0433 -3.3834 0.1442 10.0000 0.4500 -0.1493 -4.0498 0.1926 11.0000 0.5000 -0.3578 -3.9654 0.2085 12.0000 0.5500 -0.5405 -3.0262 0.1827 13.0000 0.6000 -0.6547 -1.3003 0.1142 14.0000 0.6500 -0.6656 0.9689 0.0109 15.0000 0.7000 -0.

34、5546 3.3962 0.1111 16.0000 0.7500 -0.3249 5.5195 0.2297 17.0000 0.8000 -0.0037 6.8857 0.3212 18.0000 0.8500 0.3607 7.1412 0.3644 19.0000 0.9000 0.7063 6.1110 0.3456 20.0000 0.9500 0.9675 3.8486 0.2612 21.0000 1.0000 1.0876 0.6456 0.1201表3输出解释：k为迭代次数即循环的次数，X为自变量在区间内按步长的变化，Y有两列其中第一列为计算的值第二列为计算的的值，对于P前

35、四列分别为k，X，Y最后一列为误差。我们用matlab绘制的图像如下：图54.3 四阶RungeKutta方法求解利用分量形式的程序做具体步骤如下：首先建立分量形式四阶RungeKutta函数并保存；其次建立相应的分量函数：f1 f2 并保存；建立M文件gaijineuler并保存。调用gaijineuler文件运行并输出结果：x = 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500 0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500 0.8000 0.8500 0.9000 0.9

36、500 1.0000y =0 0.0195 0.0717 0.1397 0.1998 0.2272 0.2025 0.1172 -0.0233 -0.1989 -0.3781 -0.5230 -0.5973 -0.5740 -0.4418 -0.2097 0.0928 0.4193 0.7135 0.9187 0.9887z= 0 0.7577 1.2745 1.3668 0.9530 0.0777 -1.0931 -2.3016 -3.2511 -3.6693 -3.3706 -2.3043 -0.5780 1.5495 3.7076 5.4770 6.4689 6.4036 5.1739

37、2.8805 -0.1687wucha =0 0.0195 0.0522 0.0680 0.0601 0.0274 0.0247 0.0853 0.1405 0.1757 0.1792 0.1449 0.0743 0.0234 0.1322 0.2321 0.3025 0.3265 0.2942 0.2052 0.0700图6采用向量形式的程序：输出介于K，X，Y，误差最后均在P中含有所以我们在此只给出P的值：P = 1.0000 0 0 0 0 2.0000 0.0500 0.0195 0.7577 0.0195 3.0000 0.1000 0.0717 1.2745 0.0522 4.00

38、00 0.1500 0.1397 1.3668 0.0680 5.0000 0.2000 0.1998 0.9530 0.0601 6.0000 0.2500 0.2272 0.0777 0.0274 7.0000 0.3000 0.2025 -1.0931 0.0247 8.0000 0.3500 0.1172 -2.3016 0.0853 9.0000 0.4000 -0.0233 -3.2511 0.1405 10.0000 0.4500 -0.1989 -3.6693 0.1757 11.0000 0.5000 -0.3781 -3.3706 0.1792 12.0000 0.5500

39、 -0.5230 -2.3043 0.1449 13.0000 0.6000 -0.5973 -0.5780 0.0743 14.0000 0.6500 -0.5740 1.5495 0.0234 15.0000 0.7000 -0.4418 3.7076 0.1322 16.0000 0.7500 -0.2097 5.4770 0.2321 17.0000 0.8000 0.0928 6.4689 0.3025 18.0000 0.8500 0.4193 6.4036 0.3265 19.0000 0.9000 0.7135 5.1739 0.2942 20.0000 0.9500 0.91

40、87 2.8805 0.2052 21.0000 1.0000 0.9887 -0.1687 0.0700图74.4 四阶Adams预测校正系统方法求解我们首先建立Adams预测校正系统程序dAdamsyx文件并保存；其次建立运行时的M文件Adamsd，保存被调用；输出结果为：介于输出时k，X，Y，wucha在M中均包含输出只是为了在Adams文件内画图方便所以在这我们仅给出输出量M。M =1.0000 0 0 0 0 2.0000 0.0500 0.0195 0.7577 0.0195 3.0000 0.1000 0.0717 1.2745 0.0522 4.0000 0.1500 0.13

41、97 1.3668 0.0680 5.0000 0.2000 0.1998 0.9556 0.0600 6.0000 0.2500 0.2274 0.1983 0.0277 7.0000 0.3000 0.2100 -0.7757 0.0174 8.0000 0.3500 0.1413 -1.7639 0.0687 9.0000 0.4000 0.0273 -2.5470 0.1141 10.0000 0.4500 -0.1161 -2.9248 0.1433 11.0000 0.5000 -0.2638 -2.7618 0.1478 12.0000 0.5500 -0.3871 -2.020

42、1 0.1233 13.0000 0.6000 -0.4582 -0.7759 0.0711 14.0000 0.6500 -0.4562 0.7864 0.0020 15.0000 0.7000 -0.3720 2.4011 0.0842 16.0000 0.7500 -0.2111 3.7663 0.1609 17.0000 0.8000 0.0060 4.5990 0.2171 18.0000 0.8500 0.2463 4.6903 0.2403 19.0000 0.9000 0.4691 3.9504 0.222820.0000 0.9500 0.6332 2.4351 0.1641

43、图85 结果及误差分析在Y的输出结果中第一列为为所使用的计算方法计算出的数值解，第二列为y导数值，P(M)的第一列为k，第二列为x的变化，第三列为计算的数值解y，第四列为导数值，第五列为误差值(采用分量的那个程序除外)。通过误差我们发现在某些点处的误差比较大但是均没有超过0.4。这说明改进Euler计算的结果还是比较精确的。误差分析：误差来源主要有两个：第一是方法误差例如在用RungeKutta计算式进行计算时公式本身存在的误差为，其他方法也存在类似的误差；第二个是机器误差，在没有设置字节的长度时机器自动截取小数点后四位，所以存在截断误差。在这两个误差存在导致了最终的误差。小结本次课程设计所做

44、的内容是关于微分方程数值解的问题，所涉及的方法较多在此仅选取几个较为典型的方法做了简单的阐述，由于个人能力有限所述问题或许没有使读者明白怎么回事，在此希望老师给予多多扶正。此外本文并没有就误差做详细说明尤其是方法误差希望读者自行学校解决。通过这次课程设计，我刻理解了上述的几种方法，并且对推到过程有了新的理解。在做课程设计时我初步掌握了科学论文的写作格式，和学习知识的求实务真精神。对待学术知识的求知在这次课程设计中我深有体会例如龙格库塔方法，这种创新精神是我们值得学习和推崇的。 参考文献1黄明游，冯果忱数值分析（下册）M 北京：高等教育出版社，2008：203-2072李庆扬，王能超，易大义数值

45、分析（第5版）M 北京：清华大学出版社，2008：286-2903李庆扬，王能超，易大义数值分析（第5版）M 北京：清华大学出版社，2008：297-3034 任玉杰数值分析及其MATLAB实现（MATLAB 6.X,7.X版）（附学习光盘）M 北京：高等教育出版社，2007：830-870附录：附录1改进Euler方法(向量形式的程序)程序文件：function k,X,Y,wucha,P=gjeuler(dydx,a,b,CT,h)%输入为：dydx为；为所求积分的区间；CT为初始值向量。%输出为：k表示迭代次数；X为自变量的x变化情况；Y的第一列为算值y第二列为计算值，依此类推； P的列

46、依此k,X,Y，第五列为误差的变化。n=fix(b-a)/h); %求出迭代次数 X=zeros(n+1,1); %给出向量保存自变量x的变化Y=zeros(n+1,length(CT); %根据初值条件确定Y的维数Y的第一列保存第二列保存第三%列保存依此类推下去。X=a:h:b; Y(1,:)= CT;for k=1:nf1=h*feval(dydx,X(k),Y(k,:);f1=f1;x2=X(k);y2=Y(k,:);f2=h*feval(dydx,x2,Y(k,:)+f1);f2=f2;Y(k+1,:)=Y(k,:)+(f1+f2)/2;k=k+1;end %主体循环解释同前面分析fo

47、r k=2:n+1wucha(k)=norm(Y(k)-Y(k-1); k=k+1; %利用二泛数求误差endX=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1; wucha=wucha(1:k,:);P=k,X,Y,wucha; %把各输出量保存在P中原函数文件：function dY=dydx(X,Y)dY1=Y(2);dY2=-64.*Y(1)+16.*cos(8*X);dY=dY1;dY2;gaijineuler文件：CT=0;0;h=0.05;a=0;b=1;k,X,Y,wucha,P=gjeuler(dydx,a,b,CT,h)plot(X,Y(:,1),g-,X,Y(

48、:,2),r:) %绘制变量和图像hold onx=0:0.05:1;y=sin(8.*x).*x; plot(x,y,r*) %绘制变量真实解图像保存在同一图中便于观察hold onxlabel(轴it x); ylabel(轴it y)legend(是方程解y的曲线,是解y的一阶导数,方程的精确解) %对图像中的图像进行解释附录2采用分量形式的RungeKutta做法分量形式的四阶RungeKutta函数：Runge_kuttafunction x,y,z,wucha=Runge_kutta(a,b,y0,z0,h)%输入量：为所求积分的区间，y0为的值，z0为的值。%输出量：x为自变量x

49、的变化值，y为所求的数值积分值，z为y的导数值，wucha为误差值。x=a:h:b;y(1)=y0;z(1)=z0;n=(b-a)/h+1;for i=2:n K(1,1)=f1(x(i-1),y(i-1),z(i-1); K(2,1)=f2(x(i-1),y(i-1),z(i-1); K(1,2)=f1(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); K(2,2)=f2(x(i-1)+h/2,y(i-1)+K(1,1)*h/2,z(i-1)+K(2,1)*h/2); K(1,3)=f1(x(i-1)+h/2,y(i-1)+K(1,2)*h/2

50、,z(i-1)+K(2,2)*h/2); K(2,3)=f2(x(i-1)+h/2,y(i-1)+K(1,2)*h/2,z(i-1)+K(2,2)*h/2); K(1,4)=f1(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); K(2,4)=f2(x(i-1)+h,y(i-1)+K(1,3)*h,z(i-1)+K(2,3)*h); y(i)=y(i-1)+h/6*(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4); z(i)=z(i-1)+h/6*(K(2,1)+2*K(2,2)+2*K(2,3)+K(2,4);end %主体循环同程序分析中

51、的相同for i=2:n wucha(i)=abs(y(i)-y(i-1); %求出误差end分量形式的原函数：function f1=f1(x,y,z) % f1函数f1=z;function f2=f2(x,y,z) % f2函数f2=-64*y+16*cos(8*x);分量形式的gaijineuler文件y0=0;z0=0;h=0.05;a=0;b=1;x,y,z,wucha=Runge_kutta(a,b,y0,z0,h)Y=sin(8.*x).*x;plot(x,y,b-,x,z,gp,x,Y,r*) %在同一图像画出数值解，导数值，和精确解的图像xlabel(轴it时间x); yl

52、abel(轴it y);legend(是方程解y的曲线,是解y的一阶导数,方程的精确解)附录3向量形式的函数程序RK4z：function k,X,Y,wucha,P=RK4z(dydx,a,b,CT,h)n=fix(b-a)/h);X=zeros(n+1,1); Y=zeros(n+1,length(CT); X=a:h:b; Y(1,:)= CT;for k=1:nk1=feval(dydx,X(k),Y(k,:); x2=X(k)+h/2;y2=Y(k,:)+k1*h/2;k2=feval(dydx,x2,y2); k3=feval(dydx,x2,Y(k,:)+k2*h/2); k4=

53、feval(dydx, X(k)+h,Y(k,:)+k3*h); Y(k+1,:)=Y(k,:)+h*(k1+2*k2+2*k3+k4)/6;k=k+1;endfor k=2:n+1wucha(k)=norm(Y(k)-Y(k-1); k=k+1;endX=X(1:n+1);Y=Y(1:n+1,:);k=1:n+1;wucha=wucha(1:k,:)P=k,X,Y,wucha;Runge_kutta向量函数及其M文件CT=0;0;h=0.05;a=0;b=1;k,X,Y,wucha,P=RK4z(dydx,a,b,CT,h)x=0:0.05:1;y=sin(8.*x).*x;plot(X,Y

54、(:,1),g-,X,Y(:,2),r:,x,y,r*) xlabel(时间轴it x); ylabel(轴it y)legend(是方程解y的曲线,是解y的一阶导数,方程的精确解)%程序的输入输出以及M文件与前面的大体相同在此不重复。附录4 采用adams预测校正系统做：function k,X,Y,wucha,M=dAdamsyx(dydx,a,b,CT,h)x=a; y=CT;n=fix(b-a)/h);if n0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=ifft2(Ya);figuresubpl

55、ot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 二、理想高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);Hd=ones(size(IA);r=sqrt(f1.2+f2.2);Hd(r0.2)=0;Y=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=rea

56、l(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth低通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgrid);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i

57、,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); Butterworth高通滤波器IA=imread(lena.bmp);f1,f2=freqspace(size(IA),meshgr

58、id);D=0.3;r=f1.2+f2.2;n=4;for i=1:size(IA,1) for j=1:size(IA,2) t=(D*D)/r(i,j); Hd(i,j)=1/(tn+1); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya);figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Face

59、lighting,phong); 高斯低通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);D=100/size(IA,1);r=f1.2+f2.2;Hd=ones(size(IA);for i=1:size(IA,1) for j=1:size(IA,2) t=r(i,j)/(D*D); Hd(i,j)=exp(-t); endendY=fft2(double(IA);Y=fftshift(Y);Ya=Y.*Hd;Ya=ifftshift(Ya);Ia=real(ifft2(Ya)

60、;figuresubplot(2,2,1),imshow(uint8(IA);subplot(2,2,2),imshow(uint8(Ia);figuresurf(Hd,Facecolor,interp,Edgecolor,none,Facelighting,phong); 高斯高通滤波器IA=imread(lena.bmp);IB=imread(babarra.bmp);f1,f2=freqspace(size(IA),meshgrid);%D=100/size(IA,1);D=0.3;r=f1.2+f2.2;for i=1:size(IA,1) for j=1:size(IA,2) t=r