实验二 MATLAB数值计算 常微分方程(组)的求解.docVIP

实验二? MATLAB数值计算：常微分方程(组)的求解 一、实验目的 在物理学和工程技术上，很多问题都可以用一个或一组常微分方程来描述，因此要解决相应的实际问题往往需要首先求解对应的微分方程。在大多数情况下这些微分方程通常是非线性的或者是超越方程(比如范德堡方程，波导本征值方程等)，因此往往需要使用计算机数值求解。MATLAB作为一种强大的科学计算语言，其在数值计算和数据的可视化方面具有无以伦比的优势。在解决常微分方程问题上，MATLAB就提供了多种可适用于不同场合(如刚性和非刚性问题)下的求解器(Solver)，例如ode45，ode15s，ode23，ode23s等等。本次实验将以范德堡方程的计算和地球卫星的运行轨道的仿真为例，练习使用MATLAB的常微分方程求解器，以期达到如下几个目的： 熟悉常微分方程的求解方法，了解状态方程的概念； 能熟练使用dsolve函数解析求解常微分方程； 能熟练运用ode45、ode15s求解器分别数值求解非刚性和刚性常微分方程； 学习用求解器来绘制相图的方法。 二、实验的预备知识 1．微分方程的概念 未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为微分方程。如果未知函数是一元函数，称为常微分方程（Ordinary differential equations，odes） （1） 其中t为自变量。如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为 若上式中的系数ai (t), i=1,2,…,n均与t无关，称之为常系数。 ?? 2．常微分方程的解析解 有些微分方程可直接通过积分求解，例如一阶常系数常微分方程。有些常微分方程可用一些技巧,如分离变量法,积分因子法,常数变异法,降阶法等可化为可积分的方程而求得解析解。 线性常微分方程的解满足叠加原理，从而他们的求解可归结为求一个特解和相应齐次微分方程的通解。一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。 一阶常微分方程与高阶微分方程可以互化，已知一个n阶常微分方程（显式）： （2） 设，可将上式化为n元一阶常微分方程组： 注意y1即是要求的未知函数，式(1)、(2)、(3)等价。反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。【(3)式也称状态方程，y1, y2, …,yn称为状态变量，解变量为y 】 3．微分方程的数值解法 除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程主要依靠数值解法。考虑n阶微分方程（2）式的数值求解，它等价于一阶常微分方程组(3)式。将(3)式写成通式： Matlab各种有关的求解器(Solver)都是基于(4)式求解。其中 为n个分量的列向量(Column vector)，也即状态变量； n个分量的列向量，每个元素代表常微分方程组中右边的函数表式 方程(4)式要有确定的解必须给定初值条件(t0为初始时刻)： （7） 所谓数值解法，就是在求解区间[t0 ,tf]上寻求y(t)在一系列离散节点上的近似值。称为步长，通常取为常量h。最简单的数值解法是Euler方法。 Euler法的思路极其简单：在节点处用差商近似代替(4)式中的导数 利用(4)式导出计算公式（称为Euler格式） （8） 即 (9) 因此，只要给出初始条件(7),由(8)或(9)式反复迭代，可解出各个时刻tk（k=0,1,…,m）微分方程的近似解y(tk)。 Euler方法只有一阶精度，改进方法有二阶Runge-Kutta法、四阶Runge-Kutta法、五阶Runge-Kutta-Felhberg法和先行多步法等，这些方法可用于解高阶常微分方程（组）初值问题。数值算法的主要缺点是它缺乏物理理解。 4．解常微分方程的MATLAB命令 A. 解析解 MATLAB提供了dsolve函数用于计算常微分方程的符号解，其使用格式： S = dsolve (‘方程1’, ‘方程2’,…,’初始条件1’，’初始条件2’ …,’自变量’)? 方程用字符串表示，自变量缺省值为t。导数用D表示，2阶导数用D2表示，以此类推。S返回解析解。在方程组情形，S为一个结构体数组，它的每个域存放方程组的每一个解。 例1：求下列微分方程的解析解 。 s=dsolve(’D2y=sin(2*x)-y’,’y(