6常微分方程的求解.doc

§6常微分方程的求解 一、知识背景 常微分方程在物理学，工程技术中运用非常广泛，相当重要。许多物质运动的过程用常微分方程来描述，如：质点的加速运动、谐振动、电容充放电过程及电感通电断电等过程，因此，求解常微分方程成为很多物理问题求解的一种常用方法。 而有时很难求出解析解，但可求出常微分方程的数值解以逼近解析解，以完成对摸型的研究。 求解常微分方程数值解通常有欧拉法和龙格-贝塔法。欧拉法解法的基本思想是在小区间上用差商代替导数，并通过把小区间不断地划分求极限，从而最终得到数值解。龙格-贝塔方法基本思想也是用差商代替导数，但是该方法是在小区间内运用了微分中值定理，在小区间内多取点，再取加权平均值来构造精度更高的计算式。在MATLAB系统中，主要采用龙格-贝塔法来计算常微分方程的数值解。 图 6-1 或 我们也可写成右面方程组形式： ， 这即是一阶常微分方程初值问题的一般形式。 计算指令 —— ode23，ode45 语句格式（以ode23为例）： [ t, y ] = ode23 (‘f’,tspan,y0,tol ) 语句中各符号意义如下 f：求解的常微分方程的文件名，把方程写成函数形式并存储于m文件中。方程形式为y=f(t,y)。 tspan：输入[t0,tf],分别为自变量的初始值和最终值，为单调递增（减）的积分区间。 y0：函数的初始值矢量。 Tol：误差范围，(缺省值为0.000001) [t,f]：t是输出的时间列矢量，矩阵y的每个列矢量是解的一个分量。 例1：用求数值解方法,求解常微分方程：，初始值x(0)=1。 解：定义并储存函数文件: yfun.m function exer=yfun(t,x);exer=-x.^3; 程序文件：[t,x]=ode23(yfun,[0,1],1);plot(t,x,k+,t,x,-r) * 请改指令ode45,看图形有什么区别。 Ode23，低精度，使用Runge-Kutta法的二，三阶算法。 Ode45，中等精度，使用Runge-Kutta法的四，五阶算法。 三、求解二阶以上的常微分方程的数值解 将方程在函数中迭代表达，化成一阶微分方程组，从而在命令行中编程实现。下面以单摆模型为例，说明二阶以上常微分方程的一般解法。 例2：如下图所示，在单摆模型中，以铅直方向为平衡位置，以右边为正方向建立摆角x的坐标系。 解：依据转动定律 MATLAB要求先写成一阶常微分方程组。 令：y (1) = x，y (2)=dx/dt，则： 函数文件：（保存yfun2） function exer2=yfun2(t,y);g=9.8;length=25;exer2=[y(2);-g/length*sin(y(1))]; 程序文件： ts=[0,10];y0=[0.1745,0];[t,y]=ode45(yfun2,ts,y0); subplot(1,2,1), plot(t,y),grid, subplot(2,2,2), plot(y(:,1),y(:,2)),grid * 关于函数文件的讨论： 必须严格按照格式编写，不得随意变动！ 在函数文件中，要求两个等号左边的变量名相同，在此用的是exer2。对比原方程知：exer2(1)=dy(1)/dt，exer2(2)=dy(2)/dt。 等号右边yfun2为存盘文件名，在程序文件中，必须以这名称提取。在括号中，t为时间变量，y为在一阶常微分方程组定义的量。 * 关于程序文件的讨论： ode45为解常微分方程组的指令。 调用文件，必须用函数名yfun2调用。 时间范围从0刭10，步长为默认值。 初始角位移为0.1745，初始角速度为0。 输出结果为时间变量的值t、角位移的值y(:,1),角速度的值y(:,2)（需把语句中’；’号改成’，’号才可显示）。 最后用plot指令绘制函数曲线y(1)~t（即x~t），相图y(2)~y(1)图线。 可增添语句，了解更多的东西，（请作为练习） 如，键入size(t),size(y)---可知t,y的大小，各矢量的长度，即元素个数。 键入[t,y]，可得t,y的数值解。 * 关于物理问题的一些讨论 （1）一条曲线是x~t，另一条是x~t，后者是前者的时间变化率。 （2）可以在同一图中作一条余弦曲线，（请试作一下）再和x~t比较，改变角位移初始值，从不同曲线了解x0的影响。 （3）关于相图，角位移 和角速度 是两个动力学变量，它们构成一个相平面，相平面中的每一点代表系统的一种可能的运动状态。封闭轨道相当于周期振动。 （4）再看一下阻尼摆： 例3：同样，令：y（1）= x，y（2）= dx / dt （设k=0.1） function xexer5=y