某些常微分方程数值解法与程序实现

西华大学本科毕业论文

毕业设计说明书

题 目：某些常微分方程数值解法与程序实现 学 院： 年级、专业： 学 生： 学 号： 指导教师： 完成日期：

西华大学本科毕业论文 目 录

摘 要........................................................................................................................ 4 引 言........................................................................................................................ 6 1 准备知识.................................................................................................................... 6

1.1 微分方程的解.................................................................................................. 6 1.2 常微分方程的初值问题及其解的存在唯一性.............................................. 6 1.3常微分方程的收敛性与稳定性....................................................................... 7 1.4 MATLAB开发工具简介 ................................................................................. 8 1.5 用户图形界面GUI制作 ................................................................................ 8 2 常微分方程数值解.................................................................................................... 9

2.1 Euler法 ............................................................................................................. 9 2.2 改进Euler法 ................................................................................................. 10 2.3 Runge-Kutta法 ............................................................................................... 11 3 算法及程序流程图.................................................................................................. 15

3.1Euler方法 ........................................................................................................ 15

3.1.1Euler法的算法设计 .............................................................................. 15 3.1.2Euler法的程序流程图 .......................................................................... 16 3.2 改进Euler法 ................................................................................................. 19

3.2.1 改进Euler法的算法设计 ................................................................... 19 3.2.2 改进Euler法的程序流程图 ............................................................... 20 3.3 Runge-Kutta方法 ........................................................................................... 21

3.3.1 Runge-Kutta法的算法设计 ................................................................. 21 3.3.2 Runge-Kutta法的程序流程图 ............................................................. 23

4 基于MATLAB的程序界面设计 ........................................................................... 26

4.1 用户图形界面GUI说明及详细操作过程 .................................................. 26

4.2.1 Euler向前法 ......................................................................................... 27 4.2.2 Euler向后法 ......................................................................................... 28

西华大学本科毕业论文 4.2.3 改进Euler法 ....................................................................................... 29 4.2.4 二阶Runge-Kutta方法 ....................................................................... 31 4.2.5 三阶Runge-Kutta方法 ....................................................................... 32 4.2.5 四阶Runge-Kutta方法 ....................................................................... 33 4.2.6 综合比较.............................................................................................. 35 4.2用户图形界面源程序..................................................................................... 36 5结论........................................................................................................................... 42 总结与体会.................................................................................................................. 43 谢 辞...................................................................................................................... 43 参考文献...................................................................................................................... 44

西华大学本科毕业论文 摘 要

随着社会的不断发展与进步，微分方程逐渐在越来越多的现实生活中得到了非常广泛的应用，但是，我们根据问题建立的微分方程模型中，只有非常少的模型能够得到对应的解析解，因此，通过对常微分方程数值解的剖析，我们能够更加清楚的认识和了解常微分方程。

数值分析是通过研究、开发以及分析各种数学问题来得到数值解的算法。本文主要是阐述几种常见的常微分方程数值解的方法，其中包括Euler向前法、Euler向后法、改进欧拉法、Runge-kutta法。介绍了这几种方法的基本原理及推导过程，并且用matlab实现了这些方法，包括它们算法的思想，程序流程以及具体的实现。最后，利用MATLAB中的图形用户界面GUI的功能来实现，让用户操作起来更加的得心应手。

关键词：Euler 改进欧拉 Runge-kutta MATLAB 图形用户界面

西华大学本科毕业论文 Abstract

With the continuous development and progress of society, more and more gradually in the differential equation in real life has been applied widely, but according to the differential equation model established in the model only very few can get the corresponding analytical solution, therefore, for we can know more about the differential equation, we must study the differential the number of equations solution.

The numerical analysis is through research, development and analysis of various mathematical problems to obtain the numerical solution of the algorithm. This paper is mainly for numerical solutions of ordinary differential equations of several common, including the Euler Euler forward method and backward method, the improved Euler method, Runge-kutta method is introduced. The basic principle and the derivation of these methods, and these the method is implemented by MATLAB, including their algorithm, program flow and implementation. Finally, the use of MATLAB graphical user interface GUI function to achieve, allowing users to operate more handy.

Keywords: Euler improved Euler Runge-kutta MATLAB graphical user interface

西华大学本科毕业论文 引 言

微分方程现在在许多方面已经得到了普遍应用，许多实际问题所建立的数学模型也是微分方程的定解问题。在描述系统的动态演变时，如物种的增长和蜕变、物体的运动、电路的震动瞬间、化学反应过程等都能表示为以时间t为变量的常微分方程或方程组。在实际问题中的研究中，很多时候我们很难求出微分方程的解析解或者是求解过程过于繁琐，因此我们的研究是非常的有意义。

常微分方程是微分方程中只有一个自变量的函数，本文主要是对微分方程数值解的方法进行分析，并且利用MATLAB软件将这些理论知识进行实现，既能得到微分方程的数值解，又能将这些理论知识运用到实际生活中。

本文首先对了欧拉向前算法、欧拉向后格式以及改进欧拉格式进行了分析和总结，欧拉算法相对于其他几种方法来说更易于理解，但是欧拉向前算法的弱点在于其精度不高，所以为了提高欧拉算法的精度，提出了欧拉向后算法和改进欧拉算法。其次，本文介绍了龙格-库塔法，其中包括二阶、三阶、四阶龙格-库塔法，龙格-库塔法计算过程稳定，程序也易实现，总的来说精度比欧拉算法的高。最后，本文给出了实例来证明了这几种算法的是可行的，并且得到几种算法是正确的。

1准备知识

1.1微分方程的解

一个函数如果代入微分方程后能使方程为恒等式，则称为微分方程的解。所以，对于某个区间上的所有x来说，能使微分方程Fx,f?x?,f??x?,?f?n??x??0恒成立，那么y?f?x?为常微分方程Fx,f?x?,f??x?,?f?n??x??0的解。

????1.2 常微分方程的初值问题及其解的存在唯一性

讨论常微分的初值问题：

{y'?f(x,y)y(x0)?y0西华大学本科毕业论文 (a?x?b)（1.2.1）

或记为

{dy?f(x,y)dxy(x0)?y0

(a?x?b)

只有一些特殊形式的f(x,y)，才能找到它的解析解；对于大多数常微分方程的初值问题，只能计算它的数值解。在计算中约定y(n)表示常微分方程准确解的值，

yn表示y(n)的近似值。

由微分方程理论，为保证微分方程解的存在唯一性，通常对(1.2.1)中f(x,y)对y应满足Lipschitz条件，即对任意y1、y2存在常数L?0，有

f(x,y1?f(x,y2)?Ly1?y2

（1.2.2）

1.3常微分方程的收敛性与稳定性

(一)收敛性

若对任意固定的xi?x0?ih，当 i??(同时h?0)时，有数值解

yi?y?xi?，则称该常微分方程（1.2.1）是收敛的。

常微分方程初值问题的形式还可以表示为：

yn?1?yn?h?(xn,yn,h)

若单步法有p阶精度，且增量函数?(xn,yn,h)关于y满足Lipschitz条件，即存在

L??0

?(x,y,h)??(x,y,h)?L?y?y（1.3.1）

则其数值解的整体截断误差y(xn)?yn?O(hp)。 （二）稳定性

西华大学本科毕业论文 由于计算数值解的过程中舍入误差总会存在，因此我们需要讨论其稳定性。但由于差分方法数值稳定性问题很难作一般性讨论，一般我们都只用方程

y'??y,??0 （1.3.2）

作讨论，当计算某一步yn有舍入误差时，但从此之后在计算过程中误差不会逐步扩大，我们就称这种稳定性是绝对稳定性。

1.4MATLAB开发工具简介

MATLAB是美国MathWorks公司出品的商业数学软件，它集合了数值分析、矩阵运算、信号处理和图形显示于一体，构成了一个方便、有着良好界面的用户环境，让更多不同专业的人可以方便使用。Matlab是一个非常有用的数学软件，它拥有，它不仅能够解决我们数学中遇到的各种问题，而且可以根据所求结果绘画出相应的函数图形，方便我们直接观察。无论你身处何方，从事着什么职业，只要你遇到了数学中的难题，你都可以利用Matlab，它可以为你解决你所遇到的问题，让你不用通过复杂的计算和推导来得到结果，真正实现了让用户用起来方便。

MATLAB的基本单位是矩阵。它的表达式与数学、工程计算中常用的形式十分相似，MATLAB之所以能够得到广大用户的喜爱，主要是因为它对于用户来说非常实用。由于MATLAB的许多特点，比如：MATLAB操作起来方便，使初学者易于学习和掌握；它用有非常多的工具，能够大量满足我们的需求；它能够实现程序的可视化，更用户能够更直接的使用和观察等等。总的来说，对于大部分人来说MATLAB比起其他语言更适用，特别是对于那些对语言不是特别精通的人来讲。

1.5用户图形界面GUI制作

图形用户界面是用于用户和软件相互沟通的一个界面，换句话说就是我们输入一些数据，然后调用程序，通过程序内部运算得出结果，并且把这个得到的结果显示在界面上，从而让用户更加直观、清晰的看到结果，从而方便用户的使用。

用户则通过界面上的相关内容的提示，在相对应的文本框中输入数据，然后

西华大学本科毕业论文 通过点击相对应的按钮，让计算机产生某种特定的动作或变化，比如实现计算、绘图等。本文这个设计是通过GUI来制作一个一阶常微分方程数值解法的界面。这个界面由按钮、文本、图形等组成。

其主界面由六个按钮组成，每个按钮对应着一个相应算法程序，每个算法程序可以实现一种求常微分方程的差分方法，每个按钮要求在相应的文本中输入的内容，得到对应的图形是不同的。

2常微分方程数值解

2.1 Euler法

Euler法的思想就是：在每一个小区间上，用一条切线来代替原函数曲线。从整体上看就是用一条折线来代替解（或积分）曲线，并以此来求取一系列离散结点处函数近似值。

对于常微分方程初值问题（1.2.1），在求解区间?a,b?上作等距剖分，步长

h?b?a，记xn=xn?1?h，n?1,2,...,m。用数值差商的方法，即用差商近似导数m求解常微分方程。

由数值微分向前差商公式

f'(a)?f(a?h)?f(a)h

得

y'(xn)?y(xn?h)?y(xn)y(xn?1)?y(xn)??y(xn?1)?y(xn)?hy'(xn)hh

式（1.2.1）实际上给出

y'(x)?f(x,y(x))?y'(xn)?f(xn,y(xn)) 故得

y(xn?1)?y(xn)?hf(xn,y(xn))

再由yn?y(xn),yn?1?y(xn?1)推导出向前的Euler公式：

yn?1?yn?hf(xn?yn),n?0,1,2,... （2.1.1）

西华大学本科毕业论文 同理，利用数值微分向后差商公式，可以推导出向后的Euler公式：

yn?1?yn?hf(xn?1?yn?1),n?0,1,2,...（2.1.2）

式（2.1.2）是关于yn?1的xn非线性方程，它是隐式Euler公式。 对于欧拉公式，对y(xn?1)在作Taylor展开得到：

h2y(xn?1)?y(xn?h)?y(xn)?h*y?(xn)?y??(?n),xn??n?xn?1 （2.1.3）

2!h2Euler公式（2.1.1）是由（2.1.3）展开式中截断y??(?n)所得到，因此Euler

2!公式的截断误差为O(h2)，故Euler公式是一阶方法。

2.2 改进Euler法

为了提高精度，从另一角度来对初值问题进行研究。对微分方程：

dy?f(x,y)dx

从xn到xn?1对两边进行积分，得到：y(xn?1)?y(xn)??即： y(xn?1)?y(xn)??xnxn?1xnf(x,y)dx

xn?1f(x,y)dx （2.2.1）

xn?1xn从（2.2.1）中，可以看出，只要近似地算出积分：?b?a?f(a)?f(b)?2

h?f(xn,y(xn))?f(xn?1,y(xn?1))?2 h?f(xn,y(xn))?f(xn?1,y(xn?1))?2

f(x,y)dx，就可得到

微分方程的另一种求解格式。为了减少误差，由梯形公式：

??baf(x)dx?可以得出：

xn?1xnf(x,y)dx?进一步有：

y(xn?1)?y(xn)?再用yn?1，yn分别代替y(xn?1)，y(xn)，就得到了新的求解格式：

yn?1?yn?h?f(xn,yn)?f(xn?1,yn?1)?2 （2.2.2）

西华大学本科毕业论文 公式（2.2.2）是一个隐式方程，并且计算量大。于是我们往往先用欧拉向前公式（2.1.1）求出一个初步的近似值，记作yn?1，称之为预报值，再将yn?1代替(2.2.2)右端的y(xn?1)进行计算，得到校正值y(xn?1)，于是就有：

?yn?1?yn?hf(xn,yn)(预报)?（2.2.3） ?h)?yn?1?yn??f(xn,yn)?f(xn?1,yn?1)?(校正2?方法（2.2.3）称为改进的欧拉方法。改进欧拉法的截断误差为O(h3)，欧拉方法每一步只要对f调用一次，但是因为增加了校正过程，计算量比欧拉方法增加了一倍，付出这种代价是为了提高精度。

2.3 Runge-Kutta法

通过对Euler法和改进Euler法的对比，我们可发现Euler公式（2.1.1）仅取一个点xn的斜率f(xn,yn)作为平均斜率k*的近似值，计算精度较低。而改进Euler公式（2.2.1）则是利用了两个点xn和xn?1的斜率k1?f(xn,yn)和即k*?k2?f(xn?1,yn?1)的平均值作为平均斜率k*的近似值，

k1?k2，其中k1是通2过已知的yn值，然后再利用公式（2.1.1）所求得，此时计算精度提高到了二阶。

通过上面的分析，我们可以找出改进Euler法比Euler法计算更精确的原因。那就是在确定平均斜率时，多取了一个点的斜率值。因此，受到启发，如果想办法在区间?xn,xn?1?上多找出几个点的斜率值，然后把这几个点的斜率值做加权平均，所得到的值作为k*的近似值，因此就有可能构造出更高精确度的数值计算公式，这就是Runge-Kutta法的基本思想。 （一）二阶Runge-Kutta法

在区间?xi,xi?1?上取两点分别为xi、xi???xi??h(0???1)，以及该两点处的斜率值k1和k2的加权平均(或称线性组合)来求取平均斜率k*的近似值，即

k*?c1k1?c2k2

西华大学本科毕业论文 其中点xi的斜率值k1?y?(xi)?f(xi,yi)

点xi??的斜率值k2?y?(xi??h)?f(xi??h,yi??hk1) 令：y(xi?1)?y(xi)?h(c1k1?c2k2), (2.3.1) 其中c1和c2为待定常数。

对y(xi?1)在x?xi处进行二阶Taylor展开，则有

h2y(xi?1)?y(xi)?hy?(xi)?y??(xi)2! (2.3.2)

对k2?f(xi??h,yi??hk1)在x?xi处进行一阶Taylor展开，得

k2?y?(xi)??hy??(xi) (2.3.3) 将式（2.3.3）代入式（2.3.1）中，得

y(xi?1)?y(xi)?h(c1k1?c2k2)

?y(xi)?h?c1y?(xi)?c2(y?(xi)??hy??(xi))?

?y(xi)?h(c1?c2)y?(xi)?c2?h2y??(xi)（2.3.4） 比较式（2.3.2）和式（2.3.4），得

?c1?c2?1,??1?c?.2?2 (2.3.5) ?由于方程组（2.3.5）中有三个未知量，但只有两个方程，方程的个数比未知量的个数少，因此它存在着无穷多组解。如果去c1?c2?1。则有??1，用yi?1和2yi分别代替y(xi?1)和y(xi)，此时计算公式（2.3.1）变为

h?y?y?(k1?k2)i?i?12??k1?f(xi,yi)?k?f(x.y?hk)i?1i1?2?

这一多方程就是改进Euler公式。凡是满足条件式（2.3.5）的计算格式都统

西华大学本科毕业论文 称为二阶Runge-Kutta格式。改进Euler公式就是众多的二阶Runge-Kutta公式中的一种特殊格式。

若取c1?0,c2?1，则??算公式：

1，则得到另一种形式的二阶Runge-Kutta法的计2??yi?1?yi?hk2??k1?f(xi,yi)?hh?k2?f(xi?,yi?k1)22? (2.3.6)

其中点xi?12为区间?xi,xi?1?的中点，所以公式（2.3.6）又称为中点公式。

（二）三阶、四阶Runge-Kutta法

为了进一步提高数值计算的精度，在区间?xi,xi?1?上除了xi和

xi???xi??h(0???1)两点之外，还需要再加一点xi???xi??h(????1)，以及在这三点处的斜率值k1,k2和k3的加权平均值来求平均斜率k*的近似值，此时得到一个计算公式为

yi?1?yi?h(c1k1?c2k2?c3k3) (2.3.7) 其中c1,c2和c3为待定常数。

为了求出点xi???xi??h的斜率k3?f(xi??,yi??)，必须要先找出点xi??处在位置对应的函数值yi??，我们在区间?xi,xi?1?上可以利用二阶Runge-Kutta公式，即得到

yi???yi??h(d1k1?d2k2)

其中d1和d2为待定常数。因此就有

k3?f(xi??,yi??)?f(xi??h,yi??h(d1k1?d2k2))于是得到一个计算公式为

西华大学本科毕业论文 ?yi?1?yi?h(c1k1?c2k2?c3k3)?k?f(x,y)?1ii??k2?f(xi??h,yi??hk1)??k3?f(xi??h,yi??h(d1k1?d2k2) (2.3.8)

把k2,k3在x?xi处进行二阶Taylor展开，并代入（2.3.8）中的第一个方程式，在对y(xi?1)在x?xi处进行三阶Taylor展开，并且比较两式的系数，就可以得到参数满足的方程组

???d1?d2?1??c1?c2?c3?1?1??c??c??232?1?22?c??c?3?23????cd?132?6? (2.3.9)

方程组（2.3.9）中总共有5个方程以及7个未知量，方程个数少于未知量的个数，因此解不唯一。常见的三阶Runge-Kutta公式为

h?y?y?(k1?4k2?k3)i?i?16??k1?f(xi,yi)??k?f(x?h,y?hk)ii1?222?k?f(x?h,y?hk?2hk)ii12?3 (2.3.10)

如果还需要再进一步提高计算精度到四阶，用类似上述的构造方法就可以推导出四阶Runge-Kutta法。但因为推导过程复杂，在这里就不进行推导了，四阶Runge-Kutta法也不止一个，下面给出常用的四阶方法中最常用的一个计算公式：

西华大学本科毕业论文 h?y?y?(k1?2k2?2k3?k4),i?1i?6??k1?f(xi,yi),?hh??k2?f(xi?,yi?k1),22?hh?k?f(x?,y?k2),ii?322???k4?f(xi?h,yi?hk3). (2.3.11)

3 算法及程序流程图

3.1Euler方法

1 变量说明

表1 Euler方法的符号说明

变量 相应的说明 x0 y0 x的初始条件 y的初始条件 ? a b h N 迭代精度 区间的左端点 区间的右端点 区间的步长 迭代次数 第n步x的计算结果 第n步y的计算结果 xn yn

3.1.1Euler法的算法设计 欧拉向前算法3.1

（1）输入：a,b，函数f(x,y)，步长h，初值y0 （2）赋值：N?（3）计算

b?a,n?0,x?a,y?y0 h

西华大学本科毕业论文 xn?1?xn?hyn?1?yn?hf?xn,yn?,n?0,1,2,...N 输出?x,y?

（4）如果n?N，则n?1?n,x?x0,y?y0，转3，否则停止运行 欧拉向后算法3.2

（1）输入a,b，函数f(x,y)，步长h，初值y0，精度? （2）赋值：N?（3）计算

b?a,n?0,x?a,y?y0 hxn?1?xn?hyn?1?yn?hf?xn?1,yn?1?,n?0,1,2,...N

其中yn?1用迭代法求出，输出?x,y?

（4）如果n?N且k?N，则n?1?n,x?x0,y?y0，转3，否则停止运行 3.1.2Euler法的程序流程图

图3-1是Euler向前方法的程序流程，图3-2是Euler向后方法的程序流程。

3-1

西华大学本科毕业论文 开始 输入a,b,f?x,y?,y0,h N?b?ah n?1 x1?x0?hy1?y0?hf?x0,y0? 输出x1,y1 n?N? n?n?1x0?x1,y0?y1 结束

图西华大学本科毕业论文 开始 输入a,b,f?x,y?,y0,h,a N?b?a hn?1,k?1 x1?x0?hz?y0?hf?x0,y0?k?N? z1?y0?hf(x1,z) z1?z?a? y1?z1 输出x1,y1 n?N? n?n?1 x0?x1,y0?y1结束 图3-2

西华大学本科毕业论文 3.2改进Euler法

1 符号说明

表2改进Euler方法的符号说明

变量 相应的说明 x0 y0 x的初始条件 y的初始条件 a b h N 区间的左端点 区间的右端点 区间的步长 迭代次数 第n步x的计算结果 第n步y的计算结果 xn yn 3.2.1 改进Euler法的算法设计

为了方便编程设计，我们通常把（2.2.3）改为

??yp?yn?hf?xn,yn????yq?yn?hf?xn?h,yp???yn?1?1?yp?yq??2? 算法3.3

（1）输入：区间a,b，函数f(x,y)，步长h，初值y0 （2）赋值：N?（3）计算

b?a,n?0,x?a,y?y0 hyp?yn?hf?xn,yn?

yq?yn?hf?xn?h,yp?

西华大学本科毕业论文 1?yp?yq??y,x?h?x 2 输出?x,y?

（4）如果n?N，则n?1?n,x?x0,y?y0，转3；否则停止运行 3.2.2 改进Euler法的程序流程图

图3-3是改进Euler方法的程序流程。

开始 输入a,b,f?x,y?,y0,h N?b?a hn?1 x1?x0?hyp?y0?hf?x0,y0?yp?yq2yq?y0?hf?x1,yp? y1?输出x1,y1 n?N? n?n?1 x0?x1,y0?y1结束

图3-3

西华大学本科毕业论文 3.3 Runge-Kutta方法

1 符号说明

表3Runge-Kutta方法的符号说明

变量 相应的说明 x0 y0 x的初始条件 y的初始条件 a b h 区间的左端点 区间的右端点 区间的步长 迭代次数 第n步x的计算结果 第n步y的计算结果 N xn yn

3.3.1 Runge-Kutta法的算法设计 二阶Runge-Kutta算法3.4

（1）输入：a,b，函数f(x,y)，步长h，初值y0 （2）赋值：N?（3）计算

b?a,n?0,x0?a,y?y0 hx?x0?h K1?f?x0,y0?

hh??K2?f?x0?,y0?K1?22? ?y?y0?hK2 输出?x,y?

（4）如果n?N，则n?1?n,x?x0,y?y0，转3；否则停止.

西华大学本科毕业论文 三阶Runge-Kutta算法3.5

（1）输入：a,b，函数f(x,y)，步长h，初值y0 （2）赋值：N?（3）计算

b?a,n?0,x0?a,y?y0 hx?x0?h K1?f?x0,y0?

hh??K2?f?x0?,y0?K1?22? ?hh??K3?f?x0?,y0?K2?22?? y?y0?h?K1?4K2?K3? 6 输出?x,y?

（4）如果n?N，则n?1?n,x?x0,y?y0，转3；否则停止

四阶阶Runge-Kutta算法3.6

（1）输入：a,b，函数f(x,y)，步长h，初值y0 （2）赋值：N?（3）计算

b?a,n?0,x0?a,y?y0 hx?x0?h K1?f?x0,y0?

hh??K2?f?x0?,y0?K1?22? ?hh??K3?f?x0?,y0?K2?22??

K4?f?x0?h,y0?hK3?

西华大学本科毕业论文 y?y0?h?K1?2K2?2K3?K4? 6 输出?x,y?

（4）如果n?N，则n?1?n,x?x0,y?y0，转3；否则停止 3.3.2 Runge-Kutta法的程序流程图

图3-4是二阶Runge-Kutta方法的程序流程，图3-5是三阶Runge-Kutta方法的程序流程，图3-6是四阶Runge-Kutta方法的程序流程。

开始 输入a,b,f?x,y?,y0,h N?b?a hn?1 x1?x0?hK1?f(x0,y0) hhK2?f(x0?,y0?K1)22y1?y0?hK2输出x1,y1 n?N? n?n?1 x0?x1,y0?y1结束

图3-4

图3-5

西华大学本科毕业论文 开始 输入a,b,f?x,y?,y0,h N?b?ah n?1 x1?x0?hK1?f(x0,y0)Kf(x0?hh2?2,y0?2K1) Kx0?h2,y0?h3?f(2K2)y1?y0?h6?K1?4K2?K3?输出x1,y1 n?N? n?n?1x0?x1,y0?y1 结束

西华大学本科毕业论文

开始 输入a,b,f?x,y?,y0,h N?b?a hn?1 x1?x0?hK1?f(x0,y0)hhK2?f(x0?,y0?K1)22hhK3?f(x0?,y0?K2)22K4?f(x0?h,y0?hK3)y1?y0?h?K1?2K2?2K3?K4?6 输出x1,y1 n?N? n?n?1 x0?x1,y0?y1结束

图3-6

西华大学本科毕业论文 4基于MATLAB的程序界面设计

4.1 用户图形界面GUI说明及详细操作过程

随着现在软件的更新，大部分程序都是用图形界面来实现。图形用户界面利用窗口、按钮、文本框、文字等等对象来构成一个界面，大大方便了用户的使用。根据题目的需求，在用户图形界面在设计了几个功能按钮，每一个按钮对应一种微分方程数值解法，并且在界面显示了不同解法与精确值的比较图。详细界面设计如下：

主用户图形界面由七个按钮，一个文本组成，每一个按钮又对应一个界面，其界面由两个按钮、五个需要输入数据的文本框和一个控制图形的控件组成。在执行每种不同的数值解法时，都需要输入指定类型的数据，包括常微分方程、x的取值范围、初始值以及步长，当输入数据之后，点击确认，可以求出每步迭代之后x的值和对应y的值，根据对应的数据生成对应的图形，通过图形我们可以观察到数值解与精确值的差别。

开始运行程序时，我们直接进入的是如下用户图形开始操作界面：

图4-1

西华大学本科毕业论文 在用户图形界面中，我们可以选择不同的常微分方程的数值解法，根据我们选择不同的解法，界面将跳转到相应的界面。 4.2.1 Euler向前法

当点击Euler向前这个按钮时，程序界面将跳转到如图4-2这个界面，在Euler向前这个界面中，我们看到我们需要输入一些相关的数据，如：需要求的函数，

x的取值范围，y的初始值，以及x所取的步长。

图4-2

然后，我们求y???y?x?1,y(0)?1的数值解，x的取值范围为?0,1?，步长h=0.01，将程序得到的数值解与精确解对比得出图形，结果为图4-3，当步长取值较大时，误差较大。由对比图可以看出精确解与数值解重叠效果比较好，说明其误差小，结果是正确的。

西华大学本科毕业论文

图4-3

4.2.2 Euler向后法

当点击Euler向后这个按钮时，程序界面将跳转到如图4-4这个界面，在Euler向后这个界面中，我们看到我们需要输入一些相关的数据，如：需要求的函数，

x的取值范围，y的初始值，以及x所取的步长和代数精度。

图4-4

西华大学本科毕业论文 然后，我们求y???y?x?1,y(0)?1的数值解，x的取值范围为?0,1?，步长h=0.01，代数精度为0.001，将程序得到的数值解与精确解对比得出图形，结果为图4-5，当步长取值较大时，误差较大。由对比图可以看出精确解与数值解重叠效果比较好，说明其误差小，结果是正确的。

图4-5

4.2.3 改进Euler法

当点改进Euler这个按钮时，程序界面将跳转到如图4-6这个界面，在改进Euler这个界面中，我们看到我们需要输入一些相关的数据，如：需要求的函数，

x的取值范围，y的初始值，以及x所取的步长。

西华大学本科毕业论文

图4-6

然后，我们求y???y?x?1,y(0)?1的数值解，x的取值范围为?0,1?，步长h=0.1，将程序得到的数值解与精确解对比得出图形，结果为图4-7。由对比图可以看出精确解与数值解重叠效果比较好，说明其误差小，结果是正确的。

图4-7

西华大学本科毕业论文 4.2.4 二阶Runge-Kutta方法

当点击二阶Runge-Kutta这个按钮时，程序界面将跳转到如图4-8这个界面，在二阶Runge-Kutta这个界面中，我们看到我们需要输入一些相关的数据，如：需要求的函数，x的取值范围，y的初始值，以及x所取的步长。

图4-8

然后，我们求y???y?x?1,y(0)?1的数值解，x的取值范围为?0,1?，步长h=0.1，将程序得到的数值解与精确解对比得出图形，结果为图4-9。由对比图可以看出精确解与数值解重叠效果比较好，说明其误差小，结果是正确的。

西华大学本科毕业论文

图4-9

4.2.5 三阶Runge-Kutta方法

当点击三阶Runge-Kutta这个按钮时，程序界面将跳转到如图4-10这个界面，在三阶Runge-Kutta这个界面中，我们看到我们需要输入一些相关的数据，如：需要求的函数，x的取值范围，y的初始值，以及x所取的步长。

图4-10

西华大学本科毕业论文 然后，我们求y???y?x?1,y(0)?1的数值解，x的取值范围为?0,1?，步长h=0.1，将程序得到的数值解与精确解对比得出图形，结果为图4-11。由对比图可以看出精确解与数值解重叠效果比较好，说明其误差小，结果是正确的。

图4-11

4.2.5 四阶Runge-Kutta方法

当点击四阶Runge-Kutta这个按钮时，程序界面将跳转到如图4-12这个界面，在四阶Runge-Kutta这个界面中，我们看到我们需要输入一些相关的数据，如：需要求的函数，x的取值范围，y的初始值，以及x所取的步长。

西华大学本科毕业论文

图4-12

然后，我们求y???y?x?1,y(0)?1的数值解，x的取值范围为?0,1?，步长h=0.1，将程序得到的数值解与精确解对比得出图形，结果为图4-13。由对比图可以看出精确解与数值解重叠效果比较好，说明其误差小，结果是正确的。

图4-13

西华大学本科毕业论文 4.2.6 综合比较

当点击综合比较这个按钮时，程序界面将跳转到如图4-14这个界面，在综合比较这个界面中，我们看到我们需要输入一些相关的数据，如：需要求的函数，

x的取值范围，y的初始值，以及x所取的步长。

图4-14

然后，我们求y??x2?y,y(1)?1的数值解，x的取值范围为?1,1.5?，步长h=0.1，将程序得到的几种方法的数值解与精确解的图形，结果为图4-15。

西华大学本科毕业论文

图4-15

4.2用户图形界面源程序

这几种常微分数值解在用户图形界面中实现的主要源程序代码如下： （一）Euler向前法：

function pushbutton1_Callback(hObject, eventdata, handles) fun=inline(get(handles.edit1,'String')); xmin=str2num(get(handles.edit2,'String')); xmax=str2num(get(handles.edit3,'String')); y0=str2num(get(handles.edit4,'String')); h=str2num(get(handles.edit5,'String')); x=xmin:h:xmax; N=(xmax-xmin)/h; y(1)=y0;

西华大学本科毕业论文 for n=1:N

y(n+1)=y(n)+h*feval(fun,x(n),y(n)); end

set(handles.text13,'String',num2str(x)); set(handles.text14,'String',num2str(y)); x=x';y=y';

y1=dsolve('Dy=-y+x+1','y(0)=1','x'); ezplot(y1,[xmin,xmax]); hold on; plot(x,y,'--r');

legend('解析解','数值解'); hold off;

% hObject handle to pushbutton1 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) （二）Euler向后法：

function pushbutton2_Callback(hObject, eventdata, handles) fun=inline(get(handles.edit1,'String')); xmin=str2num(get(handles.edit2,'String')); xmax=str2num(get(handles.edit3,'String')); y0=str2num(get(handles.edit4,'String')); h=str2num(get(handles.edit5,'String')); a=str2num(get(handles.edit6,'String')); x=xmin:h:xmax; N=(xmax-xmin)/h; y(1)=y0; k=1; for n=1:N

z=y(n)+h*feval(fun,x(n),y(n));

西华大学本科毕业论文 while k

z1=y(n)+h*feval(fun,x(n+1),z); if(abs(z1-z))

z=z1; end

y(n+1)=z1; end

set(handles.text9,'String',num2str(x)); set(handles.text10,'String',num2str(y)); x=x';y=y';

y1=dsolve('Dy=-y+x+1','y(0)=1','x'); ezplot(y1,[xmin,xmax]); hold on; plot(x,y,'--r');

legend('解析解','数值解'); hold off;

% hObject handle to pushbutton2 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB

% handles structure with handles and user data (see GUIDATA

（三）改进Euler法：

function pushbutton2_Callback(hObject, eventdata, handles) fun=inline(get(handles.edit1,'String')); xmin=str2num(get(handles.edit2,'String')); xmax=str2num(get(handles.edit3,'String')); y0=str2num(get(handles.edit4,'String')); h=str2num(get(handles.edit5,'String')); N=(xmax-xmin)/h;

西华大学本科毕业论文 y(1)=y0; x(1)=xmin; x(N)=0;y(N)=0; for i=1:N

x(i+1)=xmin+i*h;

y1=y(i)+h*feval(fun,x(i),y(i)); y2=y(i)+h*feval(fun,x(i+1),y1); y(i+1)=(y1+y2)/2; end

set(handles.text8,'String',num2str(x)); set(handles.text9,'String',num2str(y)); y1=dsolve('Dy=-y+x+1','y(0)=1','x'); ezplot(y1,[xmin,xmax]); hold on; plot(x,y,'--r');

legend('解析解','数值解'); hold off;

% hObject handle to pushbutton2 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) （四）二阶Runge-Kutta方法：

function pushbutton2_Callback(hObject, eventdata, handles) fun=inline(get(handles.edit1,'String')); xmin=str2num(get(handles.edit2,'String')); xmax=str2num(get(handles.edit3,'String')); y0=str2num(get(handles.edit4,'String')); h=str2num(get(handles.edit5,'String')); x=xmin:h:xmax; N=(xmax-xmin)/h;

西华大学本科毕业论文 y(1)=y0; for n=1:N

k1=feval(fun,x(n),y(n));

k2=feval(fun,x(n)+h/2,y(n)+h/2*k1); y(n+1)=y(n)+h*k2; end

set(handles.text8,'String',num2str(x)); set(handles.text9,'String',num2str(y)); x=x';y=y';

y1=dsolve('Dy=-y+x+1','y(0)=1','x'); ezplot(y1,[xmin,xmax]); hold on; plot(x,y,'--r');

legend('解析解','数值解'); hold off;

% hObject handle to pushbutton2 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) （五）三阶Runge-Kutta方法：

function pushbutton2_Callback(hObject, eventdata, handles) fun=inline(get(handles.edit1,'String')); xmin=str2num(get(handles.edit2,'String')); xmax=str2num(get(handles.edit3,'String')); y0=str2num(get(handles.edit4,'String')); h=str2num(get(handles.edit5,'String')); x=xmin:h:xmax; N=(xmax-xmin)/h; y(1)=y0; for n=1:N

西华大学本科毕业论文 k1=feval(fun,x(n),y(n));

k2=feval(fun,x(n)+h/2,y(n)+h/2*k1); k3=feval(fun,x(n)+h,y(n)-h*k1+2*h*k2); y(n+1)=y(n)+h/6*(k1+4*k2+k3); end

set(handles.text10,'String',num2str(x)); set(handles.text11,'String',num2str(y)); y1=dsolve('Dy=-y+x+1','y(0)=1','x'); ezplot(y1,[xmin,xmax]); hold on; x=x';y=y'; plot(x,y,'--r');

legend('解析解','数值解'); hold off;

% hObject handle to pushbutton2 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA) （六）四阶Runge-Kutta方法：

function pushbutton2_Callback(hObject, eventdata, handles) fun=inline(get(handles.edit1,'String')); xmin=str2num(get(handles.edit2,'String')); xmax=str2num(get(handles.edit3,'String')); y0=str2num(get(handles.edit4,'String')); h=str2num(get(handles.edit5,'String')); x=xmin:h:xmax; N=(xmax-xmin)/h; y(1)=y0; for n=1:N

k1=feval(fun,x(n),y(n));

西华大学本科毕业论文 k2=feval(fun,x(n)+h/2,y(n)+h/2*k1); k3=feval(fun,x(n)+h/2,y(n)+h/2*k2); k4=feval(fun,x(n)+h,y(n)+h*k3); y(n+1)=y(n)+h/6*(k1+2*k2+2*k3+k4); end

set(handles.text8,'String',num2str(x)); set(handles.text9,'String',num2str(y)); x=x';y=y';

y1=dsolve('Dy=-y+x+1','y(0)=1','x'); ezplot(y1,[xmin,xmax]); hold on; plot(x,y,'--r');

legend('解析解','数值解'); hold off;

% hObject handle to pushbutton2 (see GCBO)

% eventdata reserved - to be defined in a future version of MATLAB % handles structure with handles and user data (see GUIDATA)

5结论

本次的设计主要对一阶常微分方程数值解的基本思想及数学算法原理进行了描述，分别对Euler法、改进Euler法、Runge-kutta法的基本数学方法和计算公式进行了分析。本次设计还通过MATLAB对Euler法、改进Euler法、Runge-kutta法进行了程序的实现，包括算法思想、程序框图及具体的程序。更是利用MATLAB的用户图形界面GUI实现了这几种方法的具体实现，让用户拥有了一个方便、快捷、简单的用户界面。

当然，在完成这次设计的过程中，我也遇到了不少的问题，由于之前没有用过MATLAB制作过界面，所以在这次设计界面时，总是不能够把具体的代码和界面连接起来，最后，通过上网查询资料以及老师和同学的帮助，最终让这个问题得到了解决。

西华大学本科毕业论文 总结与体会

通过这次对一阶常微分方程数值解的算法设计，我学会了许多，同时也收货了许多。首先，我对常微分方程数值解的几种方法有了较多的了解，之前，我对常微分方程数值解的几种方法的认识比较片面的。在做题时，我也仅仅只是靠这些公式来得出最后的结果。至于Euler法、改进Euler法、Runge-kutta法的基本思想，数学原理及数学公式的是不了解的，更不可能通过编写程序来实现该算法了。但通过这次设计，我对Euler法、改进Euler法、Runge-kutta法有了不一样的认识，对它们各自的基本思想，数学原理及数学公式的由来有了进一步的了解。其次，通过这次程序设计，让我知道了大部分的数学问题原理都可以通过计算机编程来解决，MATLAB是一个很好解决数学问题的学习软件。数值计算往往需要完成大量繁琐公式的计算，其计算量已经达到了相当恐怖的程度，大多数情况下仅仅靠人手写和计算的能力往往不能实现，但是，利用MATLAB对这些问题进行编程后再求解，仅仅只需要几秒钟的时间就可以得到我们想要的结果，并且还可以同时用绘画出相应的图形。由此可见，MATLAB对数学问题的求解其实是非常有优势的。

谢 辞

本次课程设计我能够顺利地完成，首先，必须要感谢指导老师郑鹏社老师，感谢他热心的指导。在做毕业设计的过程中，我不明白几种方法的原理，都是老师耐心的为我讲解，直至我听懂为止，然后，当我们完成论文的时候，也是老师一遍又一遍地为我们修改，我想如果不是他的指导，这次的毕业设计我完成起来是非常困难的。其次，我还要感谢我们同组的同学们热心的帮助，每当我在编程方便遇到困难时，他们都会为我提供一些想法，帮助我共同完成我的编程，并且在空闲的时候，我们也会讨论和交流，在这其中对我也有很大的帮助。最后，在这里再次对他们表示衷心的感谢，谢谢他们的无私，因为他们的无私使我得到了无比大的收获。

西华大学本科毕业论文 参考文献

[1]Uri M.AscherLinda R. Petzold. Computer methods for ordinary differential equations and differential-algebraic equations[M],北京: 科学出版社, 2009; [2]冯康. 数值计算方法[M], 北京: 国防工业出版社,1978

[3]张韵华.奚梅成. 陈效群. 数值计算方法与算法[M],北京: 科学出版社, 2006 [4]David K. 数值分析[M], 北京: 机械工业出版社, 2003

[5]李庆扬. 关治. 白峰杉. 数值计算原理[M], 北京: 清华大学出版杜, 2000 [6] 石博强. MATLAB程序设计与实例应用[M], 北京: 中国铁道出版社,2003 [7]周恺.浅析一阶常微分方程的解法问题[J].考试周刊,2013,27:51

[8] 刘林. 一阶常微分方程初等解法研究[J]. 河套大学学报（自然科学版）, 2006, 13: 15 [9] 周义仓. 常微分方程及其应用[M] , 北京: 高等教育出版社, 1985

[10] Joseph L Z. Introduction to scientific programming[M], Berlin:Springger-Verlag Inc, 1998 [11] John H. Numercal Methods Using MATLAB Fourth Edition[M], 北京: 电子工业大学出版社, 2005

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