数学实验 常微分方程.docVIP

实验六 常微分方程的Matlab解法 一、实验目的 了解常微分方程的解析解。 了解常微分方程的数值解。 学习掌握MATLAB软件有关的命令。 二、实验内容 一根长的无弹性细线，一段固定，另一端悬挂一个质量为的小球，在重力的作用下小球处于垂直的平衡位置。若使小球偏离平衡位置一个角度，让它自由，它就会沿圆弧摆动。在不考虑空气阻力的情况下，小球会做一定周期的简谐运动。利用牛顿第二定律得到如 下的微分方程 问该微分方程是线性的还是非线性的？是否存在解析解？如果不存在解析解，能否求出其近似解？ 三、实验准备 MATLAB中主要用dsolve求符号解析解，ode45,ode23,ode15s求数值解。 s=dsolve(‘方程1’, ‘方程2’,…,’初始条件1’，’初始条件2’ …,’自变量’) 用字符串方程表示，自变量缺省值为t。导数用D表示，2阶导数用D2表示，以此类推。S返回解析解。在方程组情形，s为一个符号结构。 [tout,yout]=ode45(‘yprime’,[t0,tf],y0) 采用变步长四阶Runge-Kutta法和五阶Runge-Kutta-Felhberg法求数值解，yprime是用以表示f(t,y)的M文件名，t0表示自变量的初始值，tf表示自变量的终值，y0表示初始向量值。输出向量tout表示节点(t0,t1, …,tn)T,输出矩阵yout表示数值解，每一列对应y的一个分量。若无输出参数，则自动作出图形。 ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。ode23与ode45类似，只是精度低一些。ode12s用来求解刚性方程组，是用格式同ode45。可以用help dsolve, help ode45查阅有关这些命令的详细信息. 四、实验方法与步骤 练习1 求下列微分方程的解析解 （1） （2） （3） 方程（1）求解的MATLAB代码为： clear; s=dsolve(Dy=a*y+b) 结果为 s =-b/a+exp(a*t)*C1 方程（2）求解的MATLAB代码为： clear; s=dsolve(D2y=sin(2*x)-y,y(0)=0,Dy(0)=1,x) simplify(s) %以最简形式显示s 结果为 s =(-1/6*cos(3*x)-1/2*cos(x))*sin(x)+(-1/2*sin(x)+1/6*sin(3*x))*cos(x)+5/3*sin(x) ans =-2/3*sin(x)*cos(x)+5/3*sin(x) 方程（3）求解的MATLAB代码为： clear; s=dsolve(Df=f+g,Dg=g-f,f(0)=1,g(0)=1) simplify(s.f) %s是一个结构 simplify(s.g) 结果为 ans =exp(t)*cos(t)+exp(t)*sin(t) ans =-exp(t)*sin(t)+exp(t)*cos(t) 练习２ 求解微分方程 先求解析解，再求数值解，并进行比较。由 clear; s=dsolve(Dy=-y+t+1,y(0)=1,t) simplify(s) 可得解析解为。下面再求其数值解，先编写M文件fun8.m %M函数fun8.m function f=fun8(t,y) f=-y+t+1; 再用命令 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) 结果见图6.7.1 图6.7.1 解析解与数值解 由图7.1可见，解析解和数值解吻合得很好。 下面我们讨论实验引例中的单摆问题. 练习3 求方程 的数值解.不妨取.则上面方程可化为 先看看有没有解析解.运行MATLAB代码 clear; s=dsolve(D2y=9.8*sin(y),y(0)=15,Dy(0)=0,t) simplify(s) 知原方程没有解析解.下面求数值解.令可将原方程化为如下方程组 建立M文件fun9.m如下 %M文件fun9.m function f=fun9(t,y) f=[y(2), 9.8*sin(y(1))]; %f向量必须为一列向量 运行MATLAB代码 clear; close; [t,y]=ode45(fun9,[0,10],[15,0]); plot(t,y(:,1)); %画随时间变化图,y(:2)则表示的值 xlabel(t),yla