2013《常微分方程课程设计》指导书 1-2.doc

第1章 引 言1.1 课程设计的意义高等学校的实践教学一般包括课程实验、综合性设计(课程设计)、课外科技活动、社会实践、毕业设计等，基本上可以分为三个层次：第一， 紧扣课堂教学内容，以掌握和巩固课程教学内容为主的课程实验和综合性设计；第二， 以社会体验和科学研究体验为主的社会实践和课外科技活动；第三， 以综合应用专业知识和全面检验专业知识应用能力的毕业设计。课程实践(含课程实验和课程设计)是大学教育中最重要也最基础的实践环节，直接影响后继课程的学习以及后继实践的质量。由于课程设计是以培养学生的系统设计与分析能力为目标，通过团队式合作、研究式分析、工程化设计完成较大型系统或软件的设计题目的，因此课程设计不仅有利于学生巩固、提高和融合所学的专业课程知识，更重的是能够培养学生多方面的能力，如综合设计能力、动手能力、文献检索能力、团队合作能力、工程化能力、研究性学习能力、创新能力等。常微分方程课程设计（Curriculum Design of the Ordinary Differential Equations）是一门继数学实验和常微分方程（ODE）之后开设的实验性课程，主要是指导性的讲解方程求解的数值方法和软件编程(如MATLAB，MATHMATIC，FORTRAN等)并实现方程的解析解与数值解可视化分析的一个集中实践教学环节。其宗旨在于培养学生运用计算机分析求解方程的能力，了解通过数学模型去解决实际问题的全过程，提高常微分方程课堂教学后的理解和应用效果，同时激发和提高同学们对于具有工程背景的科学研究的热情。课程设计不仅仅是以实现相应的程序为目标，更重要的是在完成课程设计的过程中逐步培养今后遇到问题而去解决问题的能力，培养从事计算机应用开发所需要的各种能力与素质。因此，在课程设计实施中，不仅需要完成程序并进行测试，还需要撰写相应的课程设计报告。课程设计报告不仅是对课程设计的总结，也是对软件文档写作能力的初步训练。1.2 课程设计的主要内容和方法常微分方程是现代数学的一个重要分支，是人们解决各种实际问题的有效工具，它在几何、力学、物理、电子技术、自动控制、航天、生命科学、经济等领域都有着广泛的应用。这些应用也为微分方程的进一步发展提出了新的问题，促使人们对微分方程进行更深入的研究，以便适应科学技术飞速发展的需要。微分方程的首要问题是如何求一个给定方程的同解或特解。到目前为止，人们已经对许多微分方程得出了求解的初等积分方法。例如一阶常系数微分方程可化为，两边积分可得通解为，其中为任意常数。有些常微分方程可用一些技巧，如分离变量法，积分因子法，常数变异法，降阶法等可化为可积分的方程而求得解析解(显式解)。另外，线性常微分方程的解满足叠加原理，从而其求解可归结为求一个特解和相应齐次微分方程的通解。一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。一阶常微分方程与高阶微分方程可以互化，所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。如给一个阶方程设，可将上式化为一阶方程组反过来，在许多情况下，一阶微分方程组也可化为高阶方程。但是我们必须认识到，可以求得解析解的常微分方程只是很少的一部分，除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程是无法得到显示的解析解的。比如简单的方程的解无法用初等函数来表示，而Bessel方程和Riccati方程一般就无法用初等积分法求出通解。因此在实际应用中，人们开始研究某种定解条件下的特解，如研究初值问题的解的性态，引入特殊函数表达解等。同时，随着计算机技术的发展，用近似方法来求微分方程的特解而成为应用中的主流解法，即依靠一些数值解法来处理复杂多样的常微分方程解的问题。比如考虑一阶常微分方程初值问题 (1.2.1)这也是一个科学与工程计算经常遇到的问题，由于只有很特殊的方程能用解析方法求解，而用计算机求解常微分方程的初值问题都要采用数值方法。通常我们假定(1.1)中对满足Lipschitz条件，即存在常数，使对有 (1.2.2)则初值问题(1.1)的解存在且唯一。假定(1.1)的精确解为，求它的数值解就是要在区间上的一组离散点上求的近似值。通常取，称为步长，求(1.1)的数值解是按节点的顺序逐步推进求得。首先，要对方程做离散逼近，求出数值解的公式，再研究公式的局部截断误差，计算稳定性以及数值解的收敛性与整体误差等问题。最简单的数值解法是Euler法。Euler方法只有一阶精度，改进方法有二阶Runge-Kutta法、四阶Runge-Kutta法、五阶Runge-Kutta-Felhberg法和线性多步法等等，这些方法可用于解高阶常微分方程（组）初值问题。边值问题采用不同方法，如差分法、有限元法等。当然，数值算法也存在一个主要缺点，那就是数值解缺乏实际物理理解。综述起来，常微分方程课程设计的主要内容就是针对实际中提炼出的常微分方程(组)初值问题，运用计算机实现各类方程(组)的数值求解。其讨论的主要方法有：一是常微分方程的单步法，主要包括Euler法、Runge-Kutta方法，如普通Euler法，后退Euler法、梯形法以及改进的Euler法，各级Runge-Kutta方法的导出和应用分析，比较各自优点和适应范围。二是线性多步法，其中包括线性多步法的一般公式，Adams显式与隐式方法，Adams预测-校正方法，Milne方法与 Hamming方法等等。最后，讨论一阶方程组与高阶方程数值方法。1.3 课程设计报告撰写基本要求课程设计报告书写要规范，报告的开头应给出题目、专业、班级、姓名、学号、完成日期，并包括以下内容：1问题背景以及问题要求陈述清楚陈述练习题的问题背景，问题所要求做的程序设计的任务，明确问题将要达到的目的或结果；2问题分析结合所学知识能具体分析问题的重点，所涉及的知识点或性质要有概括性的分析；3程序书写要求程序能顺利实现并解决问题，加上必要的注释，MATLAB语言优化，增强程序的可读性，可编多个程序。对于有些习题的完成没有具体的程序列出，也请按操作顺序予以描述。4调试和结果对程序或操作能给出调试分析过程，最终给出程序结果，该结果应该就是问题的解决结果。5总结或拓展简单总结解决问题所遇到的难点以及解决难点的方法，指出适当拓展的可能性，也可补充一些新观点。6参考文献 报告书中可能遇到查询资料来完善内容，所涉及到的文献资料必须在文本最后以统一格式编写出来。这也是学做科研的一种必要的训练。1.4 课程设计的组织与评分方法同学们通过了解并掌握本课程设计的基础知识，并在熟悉软件环境后，通过编程实现六大块的典型模型。每一块学完之后要求同学们自选题目，独立完成编写代码，测试，并撰写课程设计报告。课程设计的考核方法为：采用每章课程设计报告的百分制打分，成绩按六章内容总成绩取平均，最终折合为五级记分制：优、良、中、及格、不及格。每个报告的发布或上交都有明确的截止日期，对于超过截止日期的提交，将会酌情扣分。16第2章 MATLAB求解常微分方程本章重点在于介绍求解常微分方程的MATLAB命令的用法，一是常微分方程(组)的解析解的命令编辑和分析，二是常微分方程(组)的数值解的命令介绍和应用。我们面对的读者应该是已经熟悉了MATLAB的基本用法以及微分方程的基本内容。如果读者对所提的基本知识上不够熟悉，则可以参考本书提供的参考文献，学习一些有关的知识，这是应该做的也是十分必要的。这里所介绍的程序均以MATLAB6.5或7.0的版本为准，不过，这些程序基本上都能在MATLAB6.0以上的版本中运行。2.1 常微分方程解析解求法对于常微分方程的求解问题，这是常微分方程教材的重点甚至是全部，目的就是得到方程的精确解的表达式，也称为解析解，还称为常微分方程的符号解。在MATLAB中，由函数dsolve( )来解决常微分方程（组）的求解问题，从而得到方程的解析解。其具体格式如下： r = dsolve(eq1,eq2,., cond1,cond2,., v) (2.1.1)其中，eq1,eq2,.为微分方程或微分方程组，cond1,cond2,. 是初始条件或边界条件，v是独立变量，默认的独立变量是t。值得注意的是，函数dsolve用来解符号常微分方程、方程组，如果没有初始条件，则求出通解，如果有初始条件，则求出特解。 这个函数就涵盖了我们在微分方程教材中的各类解析解法所能解决的方程的求解问题。例1求解常微分方程 的MATLAB程序为 dsolve(Dy=1/(x+y),x) % 回车即输出结果ans = -lambertw(-C1*exp(-x-1)-x-1注意，系统缺省的自变量为 t ，因此这里需把自变量写明。 其中：Y=lambertw(X)表示函数关系Y*exp(Y)=X。而该问题的解析解法为：一是把方程变形为，按一阶线性方程的解法直接得到通解；二是可以通过变量代换，再分离变量法也可得通解。这个解跟MATLAB求得的符号解是一样的。例 2 求解常微分方程的MATLAB程序为 Y2=dsolve(D2y*y-Dy2=0,x) % 回车即输出结果Y2 = 0 exp(C1*x)*C2我们看到有两个解，其中一个是常数0。自己可以调试为什么不写成下面形式 Y2=dsolve(y*D2y-Dy2=0,x)例 3 求常微分方程组 (2.1.2)的通解的MATLAB程序为 X,Y=dsolve(Dx=-5*x-y+exp(t),Dy=x+3*y+exp(2*t),t)解析解结果为X = exp(-1+15(1/2)*t)*C2+exp(-(1+15(1/2)*t)*C1+1/6*exp(2*t)+2/11*exp(t) Y =-4*exp(-1+15(1/2)*t)*C2-exp(-1+15(1/2)*t)*C2*15(1/2)-4*exp(-(1+15(1/2)*t)*C1+exp(-(1+15(1/2)*t)*C1*15(1/2)-7/6*exp(2*t)-1/11*exp(t)例 4 求常微分方程组 (2.1.3)的通解的MATLAB程序为X,Y=dsolve(Dx+2*x-Dy=10*cos(t),Dx+Dy+2*y=4*exp(-2*t),x(0)=2,y(0)=0,t) 解析解结果为X = -2*exp(-t)*sin(t)+4*cos(t)+3*sin(t)-2*exp(-2*t) Y = 2*exp(-t)*cos(t)+sin(t)-2*cos(t) 2.2 常微分方程数值解求法上一节的方程求解都是常微分方程的精确解法。但是，我们知道，有大量的常微分方程虽然从理论上讲，其解是存在的，但我们却无法求出其解析解，此时，我们需要寻求方程的数值解，在求常微分方程数值解方面，MATLAB具有丰富的函数，我们将其统称为solver，其一般格式为：T,Y=solver(odefun,tspan,y0) (2.2.1)该函数表示在区间tspan=t0,tf上，用初始条件y0求解显式常微分方程，odefun为显式常微分方程中的中的，可以编辑M函数(如下几例)，也可以用inline()命令赋值。tspan为求解区间，要获得问题在其他指定点上的解，则令（要求单调递增或者递减）。solver为命令ode45，ode23，ode113，ode15s，ode23s，ode23t，ode23tb之一，这些命令各有特点。我们列表说明如下：表2.1 数值求解的常微分方程的常用求解器正如第一章所说的，一阶常微分方程组和高阶微分方程可以互化，在求解的任何高阶常微分方程都可以用替换法化为一阶形式，下面两例展现了一个完整处理常微分方程初值问题的数值算法的求解过程。例5 利用MATLAB求初值问题 的解。分析：该初值问题可试着按下列格式求解析解，但是其结果提示是得不到解析解。 y=dsolve(D3y-3*D2y-y*Dy=0,y(0)=0,Dy(0)=1,D2y(0)=-1,x)? Error using = dsolveError, (in dsolve/IC) The implicit option is not available when giving Initial Conditions.解：（1）化方程组为标准形式，即把微分方程的高阶导数写为低阶导数的算式：设 则原方程化为下列等价的方程组： 满足初值条件：其中我们取定各向量分别为，常数，而。（2）把微分方程组编成m函数文件。 如：function dy=F(t,y)dy=y(2);y(3);3*y(3)+y(2)*y(1);注意：在函数文件里，虽然写微分方程时并不同时包含参数t和y，但第一行必须包含这两个输入变量。另外，向量dy必须为列向量，即中括号内用分号隔开的矩阵式输入格式。（3）调用一个微分方程的求解函数来求得原方程初值问题的数值近似解，并可以画出解曲线图。 T,Y=solver(F,tspan,y0);其中solver：求解函数名，执行命令时必须选定好函数求解器，这里我们选定ode45；F：包含微分方程的m文件; tspan：积分的数据范围，其格式为：t0,tfinal; y0为t0时刻的初值列向量。输出参数T和Y为列向量：T为时刻向量，Y表是不同时刻的函数值。具体实现程序如下： T,Y=ode45(F,0,2.5,0;1;-1);%调用F，并选择ode45求解 plot(T,Y(:,1),T,Y(:,2),T,Y(:,3); legend(y,dy, d2y); xlabel(time T); ylabel(solution Y);作出解函数的曲线图如下（其中蓝色曲线即为数值方法得到的解曲线）：图2.1 例5中方程的函数y和y以及y的近似曲线图例 6 利用MATLAB求常微分方程的初值问题 的解析解和数值解，并分别作出解曲线图并列出数值对比表。 解： （1）首先求出其解析解：r=dsolve(D2y*(1+x2)-2*x*Dy=0,y(0)=1,Dy(0)=3,x) r = 1+x3+3*x（2）然后通过一阶微分方程组转化而建立M函数文件function yy=mydy2(x,Y) yy=Y(2);Y(2)*2*x/(1+x2); （3）运行数值算法程序 x,YY=ode45(mydy2,0 3,1;3); （4）作出两类解曲线图 r=1+x.*x.*x+3.*x; %解析解r的输入 plot(x,r,-b,x,YY（:,1）,*r); legend(r,y) xlabel(independent variable x); ylabel(real solution r，numerical solution y);作出解函数的曲线图如图2.2，其中红色星点连线视为解析解曲线，而蓝色曲线为数值解曲线，图中可以看出数值解非常接近解析解（真解）。图2.2 例6中初值问题的解析解和数值解曲线图数值解和解析真解对比列表如下：x自变量y数值解y真解0110.016745911.0502424251.0502424250.0334918191.1005130251.1005130250.0502377291.1508399771.1508399770.0669836381.2012514571.2012514570.1419836381.4288141421.4288132130.2169836381.6611686451.6611669160.2919836381.9008449821.9008438180.3669836382.1503754962.1503751670.4419836382.4122951262.4122922130.5169836382.6891314082.6891262080.5919836382.9834121492.98340840.6669836383.2976715973.297670040.7419836383.6344470643.6344423780.8169836383.9962640873.9962566640.8919836384.3856498524.3856441480.9669836384.8051392454.8051360791.0419836385.2572700095.2572637081.1169836385.7445673045.7445582851.1919836386.2695584936.269551061.2669836386.8347783926.8347732821.3419836387.4427642937.4427562021.4169836388.0960417498.096031071.4919836388.7971385578.7971291361.5669836389.5485891899.548581651.64198363810.3529302210.352919861.71698363811.2126878611.212675021.79198363812.1303903612.130378381.86698363813.108571813.108561181.94198363814.1497680214.149754682.01698363815.256505915.256490142.09198363816.4313141216.431298782.16698363817.6767263717.676711882.24198363818.9952778618.995260672.31698363820.3894960520.389476412.39198363821.8619121.861890352.46698363823.4150530523.415033742.54198363825.0514598725.051437832.61698363826.7736584426.773633862.69198363828.5841781228.584153092.76698363830.4855519830.485526772.82523772932.0266977932.026670962.88349181933.6253689133.625340482.9417459135.2827509535.28272145337.0000305237表2.2 例6中函数的数值解和解析解对比列表例 7 一个实际数学模型的求解问题：采用ODE解算指令研究围绕地球旋转的卫星轨道。解：（1）问题的形成我们容易知道，轨道上的卫星在牛顿第二定律和万有引力定律的作用下有 (2.2.2)其中，引力常数G=6.672*10-11(N.m2/kg2) , ME=5.97*1024(kg)是地球的质量。假定卫星以初速度vy(0)=4000m/s在x(0)=- 4.2*107(m)处进入轨道。问题就是求出随时间变化的轨迹，取直角坐标，轨线切线方向设出轴，垂直方向设出y轴，于是位移的大小为，则卫星轨道问题就转化为求解函数的近似解曲线的相位图的问题。（2）构成一阶微分方程组，令，则可构建方程组（矩阵型式）如下 (2.2.3)（3）根据上式，首先建立M函数dYdt.m function Yd=fun7(t,Y) % t % Y global G ME % 全局变量xy=Y(1:2);Vxy=Y(3:4); % 向量定义r=sqrt(sum(xy.2); Yd=Vxy;-G*ME*xy/r3; % 方程组向量形式表达（4） 编程，程序运行。 global G ME % G=6.672e-11; ME=5.97e24; vy0=4000;x0=-4.2e7; t0=0;tf=60*60*24*9; %时间设定范围为9天 tspan=t0,tf; % Y0=x0;0;0;vy0; % t,YY=ode45(fun7,tspan,Y0);% X=YY(:,1); % Y=YY(:,2); % plot(X,Y,b,Linewidth,2); hold on %axis(image) % XE,YE,ZE = sphere(10); % RE=0.64e7; % XE=RE*XE;YE=RE*YE;ZE=0*ZE; % mesh(XE,YE,ZE),hold off % legend(x,y相位图,地球二维截面图) xlabel(variable x(t); ylabel(variable y(t);图2.3 例7的卫星轨道平面图显示下面我们举一个求解刚性（stiff）方程的例子。在用微分方程描述的一个变化过程中，若往往又包含着多个相互 作用但变化速度相差十分悬殊的子过程，这样一类过程就认为具有“刚性”。描述这类过程的微分方程初值问题称为“刚性问题”。例如，宇航飞行器自动控制系统一般包含两个相互作用但效应速度相差十分悬殊的子系统，一个是控制飞行器质心运动的系统，当飞行器速度较大时，质心运动惯性较大，因而相对来说变化缓慢；另一个是控制飞行器运动姿态的系统，由于惯性小，相对来说变化很快，因而整个系统就是一个刚性系统。描述这种系统的微分方程，在数学上常常称为刚性方程。由于刚性系统中存在两(多)重尺度，一个尺度比另外一个尺度大很多。所导致的麻烦就是在计算中很难兼顾两者。无论你用什么样的尺度(单一尺度)都不能很好刻画解的行为：一是快变行为。而另一是慢变行为。所有这样的方程在计算的时候，稳定性条件比较苛刻。比如下面的刚性方程 (2.2.4)系数矩阵的特征值为，其解的表达式为， (2.2.5)显然，特征值对应系统的最小时间常数，它随时间增大很快地便趋于零而无足轻重，而对应了一个大时间行为。因此无论你用什么样的尺度(单一尺度)都不能很好刻画解的行为：一个是快变行为，一个是慢变行为。 例 8 (刚性方程组求解)求下面刚性微分方程的解 (2.2.6)使用dsolve可知解析解为下面求数值解。建立M文件fun8.m如下function f=fun8(t,y)f=-0.01*y(1)-99.99*y(2), -100*y(2);运行MATLAB代码clear; close; t,y=ode45(fun8,0,10,2,1); plot(t,y);text(1,1.1,y1);text(1,0.1,y2);xlabel(t),ylabel(y)图2.4 例8中刚性方程ode45命令所得数值解图2.4给人的感觉似乎是始终大于0.5，但由的解析解可知，当时,两个分量均趋于0。下降极快，；而下降很慢，(见下表2.3和图2.5)。时间域设为0,400，程序为clear; close; t,y=ode45(fun8,0,400,2,1); tstep=length(t) %求计算总步数minh=min(diff(t) %最小步长maxh=max(diff(t) %最大步长结果为tstep =48261minh =5.0238e-004maxh =0.0102可见计算太慢，需要48261步才能到达400。一方面，由于下降太快，为了保证数值稳定性，步