实验9常微分方程要点分析.pptVIP

求解方程组（12）所建立的M文件如下： function ode2 clear; close; [t,y] ode23 @fun10,[0,100],[15/pi,0] ; plot t,y :,1 ; xlabel t ,ylabel y1 %%%%%%%%%%%%%%%%%%%%%%% function f fun10 t,y f [y 2 ,-9.8*sin y 1 -0.1*y 2 ]; 运行结果见图7.9，可以看出，当摆受到周围介质粘性阻力 时，振动是衰减的。 （2）高阶常微分方程初值问题 若令 ，则可以化成形如 7 或 8 的形式： 或 （3）常微分方程组初值问题的数值解法 我们可以看到，高阶方程（组）的矩阵形式与单个方程（1）形式一样，只不过此处是向量，前边的是数量。构造算法时，完全将前边所建立的适合 1 的算法搬到 8 中即可，只是注意将相应的数量换成向量即可。下边我们构造两种算法。 ①　Euler法 或写成数量形式 ②　四阶经典R-K法 或写成向量形式 7.3.2常微分方程（组）初值问题数值解的Matlab命令 ode45是常用的常微分方程（组）初值问题数值解的Matlab命令，其调用方式为 [XOUT,YOUT] ode45 odefun,xspan,y0 ，其中xspan [x0, xn]为自变量取值范围，y0为初值条件，odefun指方程中的函数f x,y ， 可以是数量函数（求解单个方程，这时y0是数量），也可以是向量函数（求解方程组，这时y0是向量），输出量XOUT是节点向量，当求解单个方程时，输出量YOUT是对应于节点的数值解向量，当求解方程组时，输出量YOUT是矩阵，一列对应于一个函数在节点的数值解向量。 数值解命令还有ode23, ode23s，ode15s，ode113等，调用方式与ode45类似，各自特点及适用范围读者可以用help命令查阅详细信息。 示例6 求解微分方程 ，先求解析解，再求数值解，并进行比较 。 先求解析解，相应的Matlab代码为 clear; s dsolve Dy -y+t+1,y 0 1,t ; simplify s 得解析解为y t+exp -t . 用ode45求数值解，先编写M文件 %M函数fun8.m function f fun8 t,y y -y+t+1; 再运行相应的MATLAB代码　[t,y] ode45 fun8,[0,1],1 运行结果为 t 0，0.0250， 0.0500，0.0750，0.1000，0.1250，0.1500， 0.1750，0.2000，0.2250，0.2500， 0.2750，0.3000， 0.3250，0.3500，0.3750，0.4000，0.4250，0.4500，0.4750，0.5000，0.5250，0.5500，0.5750，0.6000，0.6250，0.6500，0.6750，0.7000，0.7250，0.7500，0.7750，0.8000，0.8250，0.8500，0.8750，0.9000，0.9250，0.9500，0.9750，1.0000 y 1.0000，1.0003，1.0012，1.0027，1.0048，1.0075，1.0107，1.0145，1.0187，1.0235，1.0288, 1.0346，1.0408，1.0475，1.0547，1.0623，1.0703，1.0788，1.0876，1.0969，1.1065，1.1166，1.1269，1.1377，1.1488，1.1603，1.1720，1.1842，1.1966，1.2093，1.2224，1.2357，1.2493，1.2632，1.2774，1.2919，1.3066，1.3215，1.3367，1.3522，1.3679 在同一图形中画出解析解与数值解的图形，观察其结果，相应的Matlab代码为 clear; close; t 0:0.1:1; y t+exp -t ; plot t,y ;%画解析解的图形 hold on %保留已经画好的图形，如果下面要再画图，两个图形和并在一起 [t,y] ode45 fun8,[0,1],1 ; plot t,y,ro ;%画数值解图形，用红色小圈画 xlabel t ,ylabel y 运行结果图7.6，可见，解析解和数值解吻合的很好。 图7.6 数值解与精确解 示例7 如图7.7所示，一根长为 的无弹性细线,一端固定，另一端悬挂质量为 的小球。在重力的作用下小球处于平衡位置。若使小球偏离平衡位置一定角度 ，放开它，它就会沿圆弧摆动。在不考虑空气阻力的情况下小球会做一定周期的简谐运动。利用牛顿