《数值分析》教案16.doc

3.4 常微分方程初值问题数值解的MATLAB实现 在MATLAB中，有多个求解常微分方程数值解的函数命令，在此介绍几个常用的函数。 3.4.1 求常微分方程初值问题数值解的函数 ode23和ode45是求解常微分方程数值解最常用的两个函数，命令中的“ode”是英文常微分方程“Ordinary Differential Equation”的缩写，它们都采用龙格-库塔公式进行数值求解，23和45分别表示使用是2/3阶和4/5阶龙格-库塔公式。它们的调用格式基本相同。在此仅以ode23为例来说明函数的用法。函数ode23的调用格式为： [x,y]=ode23(Fun,Tspan,y0,options) 1. 该命令适用于一阶常微分方程组，，如遇到高阶常微分方程，必须先把它们变成一阶常微分方程组，即状态方程，方可使用。 2. 输入参数“Fun”为定义微分方程组的M-函数文件名，可以在文件名前加写@，或用英文格式单引号界定文件名。 3. 在编辑调试窗口中编写一阶常微分方程组的M-函数文件时，每个微分方程的格式必须都与一致，即等号左端为待求函数的一阶导数，右端函数的变量严格以“先自变量、后函数”的固定顺序输入，的下标表示微分方程的序数。 4. 输入参数“Tspan”规定了常微分方程的自变量取值范围，它以矩阵的形式输入，表示自变量。 5. 输入参数“”表示初始条件向量，。微分方程组中的方程个数必须等于初始条件数，这是求常微分方程特解所必须的条件。 6. 输入参数“option”表示选项参数，它可由ODESET函数设置，较为复杂。 7. 输出参数为微分方程组解函数的列表（x和y都是列矩阵），它包含向量x各节点和与对应向量y的第i个分量值(即第i个方程解)，i表示节点序列数。 表3-1 微分方程初值问题的计算函数 函数名 阶的精确度 如何使用 ode45 中等 首先尝试此函数 ode23 低 低精度容差或适度刚性问题 ode113 从低至高 高精度容差 ode15s 从低至中 如果用ode45求解很慢 ode23s 低 低精度容差的线性刚性系统 ode23t 低 问题是适度刚性的 ode23tb 低 用低精度容差求解刚性系统 8. 输出参数缺省时输出解函数的曲线，即函数及其各阶导数的曲线。求解微分方程的命令还有ode45，ode113，ode15s，ode23s，ode23t，ode23tb（见表3-1）等。刚性方程组又称为Stiff方程组，其特点是出现解的分量级差别很大，给数值计算带来很大的困难。在化学反应、电子网络和自动控制等领域中经常遇到。 3.4.2 ode23与ode45使用方法举例 【例3-25】求二阶常微分方程 的通解与特解。 解：将二阶常微分方程变换成两个一阶常微分方程组： 第一种解法：先建立M-文件。 %定义输入,输出变量和函数文件名 Function dy=myfun_1(x,y) %明确dy的维数,用微分方程组时不可缺省 dy=zeros(2,1); %dy(m)表示y的m阶导数;y(n)表示y的第n列 dy(1)=y(2); %与方程组中第二个微分方程对应 dy(2)=2*y(2)-2*y(1)+5*exp(2*x)*sin(x); 在命令窗口键入： [x y]=ode23(@myfun_1,[0,1],[-2;-3]) 回车得到： x = 0 0.0533 …… 1.0000 y = -2.0000 -3.0000 -2.1627 -3.0957 …… -1.7691 12.8927 显示的结果是自变量x和两个待求函数和的对应数据。其中用“……”代替了许多输出数据。 再键入： x=size(x),y=size(y) x = 15 1 y = 15 2 表明自变量被分为15个节点，并计算出了对应点上的函数及其一阶导数的取值。 第二种解法：使用内联函数inline。 在命令窗口键入： myfun_2=inline([y(2);2*y(2)-2*y(1)+5*exp(2*x)*sin(x)],x,y); [x y]=ode23(myfun_2,[0 1],[-2 -3]) 回车得到与第一种方法相同的结果。 在编辑窗口再键入（不写输出参量）： ode23(@myfun_1,[0,1],[-2;-3]) legend(特解函数,一阶导函数) text(0.8,-2,特解函数,FontSize,9) text(0.8,3.5,一阶导函数,FontSize,9) 运行得到图3-1所示的微分方程的图示解。 图3-1 特解函数及其导数曲线 写有输出参量[x y]时，得出微分方程的数值解，