《Matlab在常微分方程求解中的应用》由会员分享，可在线阅读，更多相关《Matlab在常微分方程求解中的应用（50页珍藏版）》请在金锄头文库上搜索。
1、Matlab在常微分方程 求解中的应用实验目的（1）学会用Matlab软件求解微分方程的初值问题（2）了解微分方程数值解思想，掌握基本的微分方程数值解方法（3）学会根据实际问题建立简单微分方程数学模型（4）了解计算机数据仿真、数据模拟的基本方法 n17世纪：初等解法n18世纪：初等解法和无穷级数方法n19世纪：解的存在性、奇点理论、定性 理论、稳定性理论l 包含一个自变量和它的未知函数以及未知函数的导数的 等式l 形成和发展与力学、天文学、物理学及其他自然科学技 术的发展互相促进和推动常微分方程定理 设函数在区域上连续，且在区域D内满足李普希兹(Lipschitz)条件，即存在 正数L，使得对
2、于R内任意两点与,恒有则初值问题(1)的解存在并且唯一。常微分方程的解析解求微分方程（组）的解析解命令:dsolve(方程1, 方程2,方程n, 初始条件, 自变量)结 果：u = tg(t-c)记号：在表达微分方程时，用字母D表示求微分，D2、D3等 表示求高阶微分。D后所跟字母为因变量，自变量可以指定 或由系统规则选定为缺省。例如：微分方程 可以表示为D2y=0.解 输入命令: y=dsolve(D2y+4*Dy+29*y=0,y(0)=0,Dy(0)=15,x)结 果 为 : y =3e-2xsin（5x）解 输入命令：x,y,z=dsolve(Dx=2*x-3*y+3*z,Dy=4*x
3、-5*y+3*z,Dz=4*x-4*y+2*z, t)；x=simple(x) % 将x化简y=simple(y)z=simple(z)结 果 为：x = (c1-c2+c3+c2e -3t-c3e-3t)e2ty = (c1e-4t+c2e-4t+c2e-3t-c3e-3t+c1-c2+c3)e2tz = (-c1e-4t+c2e-4t+c1-c2+c3)e2t 虽然说解析解是最精确的，但是实际问 题中常要求研究常微分方程的数值解常微分方程的数值解常微分方程中只有一些典型方程能求出初等解(用初 等函数表示的解) 。另外，有些初值问题虽然有初等解， 但由于形式太复杂不便于应用。因此，有必要探讨
4、常微 分方程初值问题的数值解法。以下主要介绍一阶常微分方程初值问题几种经典数值 解方法：欧拉法及改进的欧拉法。其它方法：龙格-库塔法、阿达姆斯方法；一阶微分方程组与高阶方程初值问题的数值解法；二阶常微分方程值问题的差分方法等。求解常微分方程初值问题的数值解的整体思路：（1 ）寻求准确解在一系列离散节点：上的近似值 称为问题的数值解，数值解所满足的离散方程统称为差称为步长，实用中常取定步长。分格式，建立数值解法的一些途径1、用差商代替导数若步长h较小，则有故有公式：此即欧拉法Euler法的几何意义:找到了积分的一条近似折线找到了积分的一条近似折线2、使用数值积分 对方程y=f(x,y), 两边由
5、xi到xi+1积分，并利用梯形公式，有：实际应用时，与欧拉公式结合使用：此即改进的欧拉法故有公式：3、使用泰勒公式以此方法为基础，有龙格-库塔法、线性多步法等方 法。4、数值公式的精度当一个数值公式的截断误差可表示为O（hk+1）时 （k为正整数，h为步长），称它是一个k阶公式。k越大，则数值公式的精度越高。 欧拉法是一阶公式，改进的欧拉法是二阶公式。龙格-库塔法有二阶公式和四阶公式。线性多步法有四阶阿达姆斯外插公式和内插公式。ODE 指令列表 nMATLAB 用於求解常微分方程式的指令: 指令方法应用ODE类别ode45Explicit Runge-Kutta (4, 5) pair of
6、Dormand-PrinceNonstiff ODEode23Explicit Runge-Kutta (2, 3) pair of Bogacki and ShampineNonstiff ODEode113Variable order Adams-Bashforth-Moulton PECE solverNonstiff ODEode15sNumerical differentiation formulas (NDFS)Stiff ODEode23sModified Rosenbrock formula of order 2Stiff ODEode23t Trapezoidal rule
7、with a “free” interpolantStiff ODEode23tbImplicit Runge-Kutta formula with a backward differentiation formula of order twoStiff ODE适用于 Nonstiff 系統 n速率（即微分值）差异相常大 n使用一般的ode45、ode23或ode113來求解,可能会使 得积分的步长(Step Sizes)变得很小,以便降低积分 误差至可容忍范围以內,会导致计算时间过长 n专门对付Stiff系统的指令，例如ode15s、ode23s、 ode23t及ode23tb 指令指令项目
8、繁多，项目繁多，最主要可分最主要可分两大类：两大类：适用于 Nonstiff 系統 一般的常微分方程式都是 Nonstiff 系統 直接采用 ode45、ode23 或 ode113 來求解 提示 n使用 Simulink 來求解常微分方程式Simulink是和MATLAB共同使用的一套软件 可使用拖拉的方式來建立动力系统 可直接产生C语言代码或进行动画演示 功能非常强大 ODE ODE 指令基本用法指令基本用法 n使用ODE指令时，必须先将要求解的ODE表示成一个函数 输入为t(时间)及y(状态函数: State Variables） 输出则为dy（状态变量的微分值）ODE ODE 指令基本
9、用法指令基本用法 n若ODE函数的文件为odeFile.m,则调用ODE指令的格式如下： t,y = solver(odeFile,t0,t1,y0) t0,t1 是积分的时间区间 y0代表起始条件（Initial Conditions） solver是前表所列的各种ODE指令 t是输出的时间向量 y是对应的状态变量向量 t,y=solver（f,ts,y0,options）ode45 ode23 ode113 ode15s ode23s由待解 方程写 成的m- 文件名ts=t0 , tf， t0、tf为自 变量的初 值和终值函数的 初值ode23：组合的2/3阶龙格-库塔-芬尔格算法 ode
10、45：运用组合的4/5阶龙格-库塔-芬尔格算法自变 量值函数 值用于设定误差限(缺省时设定相对误差10-3, 绝对误差10-6), 命令为：options=odeset（reltol,rt,abstol,at）, rt，at：分别为设定的相对误差和绝对误差.1、在解n个未知函数的方程组时，y0和y均为n维向量， m-文件中的待解方程组应以y的分量形式写成.2、使用Matlab软件求数值解时，高阶微分方程必须 等价地变换成一阶微分方程组.注意:Van der Pol微分方程n1928年荷兰的范德波耳(Van der Pol)为描述LC 回路的电子管振荡器建立了著名的vanderPol方 程.它在
11、自激振荡理论中有着重要的意义,一直作 为数学物理方程中的一个基本方程.这是一个具 有可变非线性阻尼的微分方程,代表了一类极为 典型的非线性问题.和其他非线性微分方程在数 学上无法精确求解一样,人们一直在努力寻找求 解这类方程近似解析解的方法,并乐于用Van der Pol方程来检验求解方法的有效性. Van der Pol微分方程其方程形式为：n令将Van der Pol微分方程化成标准形式 写成向量的形式：为一个向量，代表状态变量n设 =1，ODE文件(vdp1.m)显示如下： type vdp1.m function dy = vdp1(t, y)mu = 1;dy = y(2); mu*
12、(1-y(1)2)*y(2)-y(1); 建立了vdp1.m后，即可选用前述ODE指令来求解 n在 =1时，van der Pol方程并非 Stiff系统， 所以使用ode45来显示结果 odeBasic01.mode45(vdp1, 0 25, 3 3); 0 25代表积分的时间区间，3 3 则代表 起始条件 因为没有输出变量，所以上述程序执行结束后， MATLAB 只画出状态变量对时间的图形 odeGetData01.mt, y = ode45(vdp1, 0 25, 3 3);plot(t, y(:,1), t, y(:,2), :);xlabel(Time t); ylabel(Sol
13、ution y(t) and y(t);legend(y(t), y(t);n也可以画出 及 在相位平面（Phase Plane）的运动情况 odePhastPlot01.mt, y = ode45(vdp1, 0 25, 3 3);plot(y(:,1), y(:,2), -o);xlabel(y(t); ylabel(y(t);n当 值越来越大时，Van der Pol方程就渐 变为一个Stiff系统，此时若要求解该类问 题，就必须用对应于Stiff系统的指令 n将 值改成 1000， ODE 文件改写成（ vdp2.m）： type vdp2.m function dy = vdp2(t
14、, y)mu = 1000;dy = y(2); mu*(1-y(1)2)*y(2)-y(1); n选用专门用来解决 Stiff 系统问题的 ODE 指令，例如 ode15s，來求解此系统并作 图显示 ode15s01.mt, y= ode15s(vdp2, 0 3000, 2 1); subplot(2,1,1); plot(t, y(:,1), -o); xlabel(Time t); ylabel(y(t); subplot(2,1,2); plot(t, y(:,2), -o); xlabel(Time t); ylabel(y(t);n由上图可知， 的变化相当剧烈（超过 ），这也就是
15、 Stiff 系统的特点 n绘制二维平面相位图：ode15s02.mn若要产生在某些特定时间点的状态变量值， 则调用 ODE 指令的格式可改成：t,y = solver(odeFile, t0,t1, ,tn, y0)其中t0, t1, , tn即是特定时间点所形成的向量 t, y= ode15s(vdp2, 0 3000, 2 1);subplot(1,1,1);plot(y(:, 1), y(:, 2), -o);xlabel(y(t); ylabel(y(t)ODE 指令的选项nODE 指令可以接受第四个输入变量，代表积分过程用到的各种选 项（Options），此种ODE指令的格式为： t,y = solver(odeFile, t0,tn, y0, options)其中 options 是由 odeset 指令來控制的结构变量 结构变量即包含了积分过程用到的各种选项 odeset的一般格式如下： options = odeset(name1, value1, name2, value2, ) nname1的值为value1，name2的值为value2 n也可以只改变一个现有的options中某栏
