* 2010-05-10 收到, 2010-08-13 改回
* * 金晓龙, 男, 1965年生, 硕士,副教授, 研究方向:应用电子技术。
文章编号: 1003-5850( 2010) 10-0026-03
Excel 求解常微分方程组
Excel Solution of Ordinary Differential Equations
金晓龙
(广东女子职业技术学院 广州 511450)
【摘 要】在建立实际问题的数学模型时, 常需要建立各物理量随时间变化的常微分方程组及求常微分方程组数
值解, 介绍借助 Excel的工作表和自定义宏函数, 用经典的四阶龙格-库塔法求常微分方程组初值问题的数值解。
【关键词】Excel, 常微分方程组, 数值解, 龙格-库塔法
中图分类号: TP311 文献标识码: A
ABSTRACT In establishing t he the mathematical m odel o f practica l pr oblems, it o fen need t o establish the o rdinar y differ ential
equations w hich phy sical quantities change over time , so lve the ordinar y differ ent ial equations fo r the numerical solutions. This
ar ticle describes how t o solve the ordinar y differ ential equat ions fo r the numer ical so lutions of initial v alues problems, with Excel
w orksheet and custom macro functions, using the classica l Runge-Kut ta method.
KEY
word
word文档格式规范word作业纸小票打印word模板word简历模板免费word简历
S Excel, o rdinar y differ ential equations, numer ical so lution, Runge-Kutta m ethod
在科学应用中, 常常需要建立实际问题的数学模
型,建立各种变量的常微分方程组及求解常微分方程
组。微分方程是指方程中未知的是一个变量或几个变
量的函数,并且在方程中不仅有函数本身而且有它们
的导数。如果方程中的未知函数是有多个变量的函数,
该方程为偏微分方程;若未知函数是只有一个变量的
函数,该方程为常微分方程。微分方程中出现的未知函
数导数的最高阶数,叫做微分方程的阶数。在实际应用
中,经常需要建立若干个物理量随时间变化的常微分
方程组,求常微分方程组初值问题的数值解是经常遇
到的问题, 本文介绍使用 Excel及龙格-库塔法对该问
题求解。
Excel是微软的电子表格软件, 广泛应用于经济、
金融、
工程
路基工程安全技术交底工程项目施工成本控制工程量增项单年度零星工程技术标正投影法基本原理
、科研、教学等领域,各单位通常使用该系统
从事日常的数据处理、绘制图表等工作。正因为 Excel
应用广泛,而求解常微分方程组的数值解问题经常遇
到,所以,在 Excel上对常微分方程组求解问题具有一
定的实际意义。Excel中的自定义宏函数是一段定义
好的操作,可以是一段程序代码,也可以是一连串的指
令集合。宏的作用是可以自动执行重复操作,从而节省
时间,提高工作效率。使用 Excel中的自定义宏函数,
在其中编写龙格-库塔具体算法,然后在 Excel表格中
调用该函数, 求解微分方程组,将步长放置在 Excel单
元格中,需要改变步长时,只需要修改相应
单元
初级会计实务单元训练题天津单元检测卷六年级下册数学单元教学设计框架单元教学设计的基本步骤主题单元教学设计
格中数
值,利用 Excel的图表功能画出数值解的模拟曲线。
1 常微分方程组求解
在工程实际与科学研究中遇到的微分方程往往比
较复杂,在很多情况下, 都不能给出解析表达式, 这些
情况下不适宜采用解析法来求解,而需采用数值解法
来求近似解。常微分方程数值解法的思路是:对求解区
间进行剖分, 然后把微分方程离散成在节点上的近似
公式
小学单位换算公式大全免费下载公式下载行测公式大全下载excel公式下载逻辑回归公式下载
或近似方程, 最后结合定解条件求出近似解。
常微分方程组与常微分方程解法不同之处在于用
向量形式构成龙格-库塔公式。对一阶常微分方程组的
初值问题。
y
′
i= f i ( x , y 1 , y 2,⋯, y N ) , ( i= 1, 2, ⋯, N )
y i ( x 0) = y 0i , ( i= 1, 2, ⋯, N )
若采用向量的记号: y = ( y1 , y 2 , ⋯, y N ) T , y 0 =
( y 01 , y 02,⋯, y 0N ) T , f = ( f 1 , f 2 ,⋯, f N ) T。
则上述方程组的初值问题可表示为:
y′= f ( x , y )
y ( x 0 ) = y 0
求解这一初值问题的四阶龙格-库塔公式为:
y n+ 1= y n+
h
6 ( k
1 , 2k2 , 2k3+ k4 )
k1= f ( x n, y n )
k2= f ( x n+
h
2
, y n+
h
2
k1 )
k3= f ( x n+
h
2
, y n+
h
2
k2 )
k4= f ( x n+ h, y n+ hk 3)
·26· (总 786) Excel求解常微分方程组 2010 年
高阶常微分方程(或方程组)可以化为一阶常微分
方程组,对下列 m 阶微分方程
y
( m )= f ( x , y , y′,⋯, y( m- 1) )
初始条件为
y ( x 0 ) = y 0 , y′( x 0 ) = y′0 ,⋯, y( m- 1) ( x 0) = ym- 10
引进新的变量
y 1= y , y 2, = y′,⋯, y m= y ( m- 1)
可将m 阶方程化为一阶方程组
y
′
1= y 2
y
′
2= y 3
⋯
y
′
m- 1= y m
y
′
m= f ( x , y 1 , y 2,⋯, y m)
初始条件相应转化为
y 1( x 0) = y 0, y 2 ( x 2) = y
′
0, ⋯, y m( x 0 ) = y ( m- 1)0
2 在 Excel 中使用龙格-库塔公式
对常微分方程组
y
′
1= xy 1+ y 2
y
′
2= xy 2
y 1 ( 0) = 1
y 2 ( 0) = 1
在区间 0≤x≤1求
解,该方程组解析解为 y 1= ( x+ 1) e
1
2
x 2
y 2= e
1
2 x
2
,在 Excel中使
用龙格-库塔公式求该方程组的数值解, 比较数值解与
解析解之间的误差。
编写宏函数。启动Excel, “工具”→“宏”→“Visual
Basic编辑器”, 在资源管理器窗口, 如图 1右上角, 右
击某一图表, 选“插入”→“模块”,在该模块填入图 1所
示代码。注意: 当创建将用在工作表公式中的自定义
函数时,要确保这些代码位于普通的 VBA 模块中, 如
果将自定义函数放在 Sheet 或者 T hisWorkbook 的代
码模块中,那么在公式中就不能运行。图 1中代码实现
龙格-库塔公式运算,由于是方程组,需要求两个解 y1、
y 2, 所以在调用函数时使用了一个逻辑变量, 为 Tr ue
返回 y 1 的值, 为 False返回 y 2的值。
代码中的 k n、l n为向量 kn中成员。
在 Excel单元格调用宏函数。由于要用前一个状
态的 x、y 1、y 2、步长 h来计算下一个状态的 y1、y2 , 在
A6单元格输入= ( A5+ $B$1) ,在 B6单元格输入=
rr($B$1, A5, B5, C5, True) , 将步长 h 放在 B1单元
格有利于调节步长重新计算, 在 C6单元格输入= =
rr($B$1, A 5, B5, C5, False) ,其余单元格用 A6、B6、
C6填充,在式子中$B$1使用绝对地址、其余使用相
对地址是为了完成填充。当下次再打开该 Excel文件
图 1 实现龙格-库塔法的函数
时,可能会禁止使用宏, 可以采用两种解决方法, 一是
设置宏安全性, 选“工具”→“宏”→“安全性”→“中”,
关闭并重新打开文件,当提示允许未签署的宏运行时
单击“启用”, “安全性”为高”或“非常高”宏未经签署,
被应用程序自动禁用,如果您已验证可以信任未签署
的宏的源,可以将安全性设为“中”或“低”;二是使用
M icrosof t Authent icode 技术允许通过使用数字证书
对文件或宏工程进行数字签名,用来创建此签名的证
书将确认宏或文档是否来源于签名者, 并且签名确认
它未被改变, 当设置宏的安全级别时,可根据宏是否由
可靠来源列表上的开发者进行数字签名来运行宏。
在单元格填充输入公式如图 2, 经运算后的结果
如图 3,可以看到数值解与解析解非常接近。
图 2 填充公式
图 3 运行结果
利用Excel的绘图功能,显示数值解的模拟曲线。
选“插入”→“图表”→“XY 散点图”→“平滑线散点
图”,下一步,“数据区域”选 A5: E15,选“系列产生在
列”, 产生 4 个序列, 分别命名“y 1 数值解”、“y 2数值
解”、“y 1解析解”、“y 2解析解”,其中 y 1解析解、y 2解
(下转第 30 页)
·27· 第 23 卷 第 10 期 电 脑 开 发 与 应 用 (总 787)
正确的,就可以继续操作, 如果三次都错了,通常卡就
会被锁定。那么,这个密码校验程序又是如何编程实现
呢? 利用生活中熟悉的问题有意思的引导学生使用多
种方法实现。
程序
设计
领导形象设计圆作业设计ao工艺污水处理厂设计附属工程施工组织设计清扫机器人结构设计
学习到这个阶段,一定要给学生建立算
法的概念。所谓算法, 为解决问题采取的方法和步
骤[ 1]。不管学生拿到一个什么样的问题,首先必须在头
脑里有非常清晰的解决这个问题的思想和步骤。这个
与任何程序设计语言无关, 与学生已有的认知结构密
切相关。学生可以用自然语言描述算法,也可以用流程
图来描述。相比自然语言,使用流程图描述算法结构更
清晰。所谓编程,无非就是用编译器能够理解的语言把
解决问题的具体步骤即流程图“翻译”出来,让计算机
能够理解和执行。所以,算法是非常重要的,在编写程
序前,应该先逻辑清晰地描述解决这个问题的算法, 然
后再用相应的语言环境编写程序,这样有助于初学者
理解程序设计,尤其是循环结构的程序设计。
1. 3 从算法分析的高度总结
前边谈到,在程序设计教学中一定要突出解决问
题算法的描述。不光在编写程序时要做到心中有明晰
的解决某一个具体问题的算法流程图; 同时还要强调
算法分析与总结,即进行某一类问题的算法分析而不
是一个问题本身的分析。每一类问题的解决,都至少有
一种算法能解决它, 在教学中要从算法的高度来进行
总结。正如“授人以鱼,不如授人以渔”,教学的目的就
是要使学生能把习得的内容迁移到新情境中去。知识
越具体,应用的范围就越狭窄,只能用于个别的情境,
也容易遗忘; 概括性越高, 其应用范围就越广,随时可
用于任何情境中类似的问题,也利于保持, 有助于学习
的迁移。教学中将循环结构常见问题如累加、连乘、穷
举、迭代、递归、求最大值或最小值、查找、排序、回溯以
及作图题(如输出杨辉三角形、输出九九乘法表)和有
关矩阵的问题(如求鞍点)等从算法的高度总结归纳,
有利于帮助学生建构思想方法层次上的编程观念。
2 总 结
程序设计是一门比较枯燥的课程, 如何让学生学
得轻松有趣是教师们值得研究的问题。所以笔者基于
建构主义和情境创设的教学理念,创设合适的问题情
境导入新课, 结合生活中学生熟悉或感兴趣的实例进
行课程主题讲述, 最后从算法高度进行总结,帮助学生
建构思想方法层次上的编程观念。实践教学证明,该教
学方法易于学生接受和掌握。通过这种教学方法进行
教学,取得了较好的教学效果,学生对循环结构的理解
更加透彻,应用更加熟练。
参 考 文 献
[ 1] 谭浩强. C 语言设计(第 3 版) [ M ] . 北京: 清华大学
出版社, 2005.
[ 2] Herber t Sch ildt. C 语言大全(第 2 版) [ M ] . 北京: 电
子工业出版社, 1994.
[ 3] 何克抗.建构主义——革新传统教学的理论基础[ J] .
电化教育研究, 1997, 17( 3) : 3-9.
[ 4] 高 楠. 艺术心理学(第 1 版) [ M ] . 沈阳: 辽宁人民出
版社, 2004.
[ 5] 黄玲玲, 阳小华.程序设计教学中的情境创设[ J] . 计算
机教育, 2007( 8) : 34-36.
[ 6] 胡正国, 吴 健, 邓正宏. 程序设计方法学 [ M ] . 北
京: 国防工业出版社, 2003.
(上接第 27 页)
析解是连续变化量, y 1数值解、y 2数值解是离散量, 且
前者大于后者,下一步, 图表标题为“常微分方程组数
值解图示”,完成, 如图 4所示。
图 4 数值解与解析解的比较图示
3 结束语
建立及解常微分方程组是实际应用中常见问题,
有些微分方程组难以求或没有解析解,通常求数值解,
这是应用计算机解决问题的常见方法, 虽然可以通过
Matlab 或其他语言求解微分方程组,由于数据通常放
在 Excel 中, 使用 Excel求解更方便、灵活,而且通过
编程可以更好地认识求解的过程。
参 考 文 献
[ 1] 樊自安. 二元常系数微分方程组解的表达式[ J] . 石家
庄学院学报, 2009, 10( 5) : 35-38.
[ 2] 程开敏. 常线性微分方程组的求解方法[ J] . 重庆三峡
学院学报, 2008, 23( 3) : 76-79.
[ 3] 吴幼明. 一类常系数微分方程组的通解[ J] . 四川理工
学院学报, 2006, 20( 4) : 68-71.
[ 4] 曹玉平. 一阶线性常系数微分方程组的矩阵解法[ J] .
河北理工学院学报, 2004, 23( 2) : 104-107.
[ 5] 吴幼明. 一类二阶常系数微分方程组的通解[ J] . 佛山
科学技术学院学报, 2002, 19( 6) : 10-14.
·30· (总 790) 循环结构程序设计教学方法 2010 年