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