数学实验-常微分方程VIP专享VIP免费

数学实验-常微分方程1/9实验六常微分方程的Matlab解法一、实验目的1．了解常微分方程的解析解。2．了解常微分方程的数值解。3．学习掌握MATLAB软件有关的命令。二、实验内容一根长l的无弹性细线，一段固定，另一端悬挂一个质量为m的小球，在重力的作用下小球处于垂直的平衡位置。若使小球偏离平衡位置一个角度，让它自由，它就会沿圆弧摆动。在不考虑空气阻力的情况下，小球会做一定周期的简谐运动。利用牛顿第二定律得到如下的微分方程0)0(',)0(,sin"0mgml问该微分方程是线性的还是非线性的？是否存在解析解？如果不存在解析解，能否求出其近似解？三、实验准备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。可以用helpdsolve,helpode45查阅有关这些命令的详细信息.四、实验方法与步骤练习1求下列微分方程的解析解（1）bayy'（2）1)0(',0)0(,)2sin(''yyyxy（3）1)0(',1)0(',','gffgggff方程（1）求解的MATLAB代码为：clear;s=dsolve('Dy=a*y+b')结果为数学实验-常微分方程2/9s=-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)练习２求解微分方程,1)0(,1'ytyy先求解析解，再求数值解，并进行比较。由clear;s=dsolve('Dy=-y+t+1','y(0)=1','t')simplify(s)可得解析解为tety。下面再求其数值解，先编写M文件fun8.m%M函数fun8.mfunctionf=fun8(t,y)f=-y+t+1;再用命令clear;close;t=0:0.1:1;y=t+exp(-t);plot(t,y);%化解析解的图形holdon;%保留已经画好的图形，如果下面再画图，两个图形和并在一起[t,y]=ode45('fun8',[0,1],1);plot(t,y,'ro');%画数值解图形，用红色小圈画xlabel('t'),ylabel('y')结果见图6.7.100.20.40.60.8111.051.11.151.21.251.31.351.4ty图6.7.1解析解与数值解数学实验-常微分方程3/9由图7.1可见，解析解和数值解吻合得很好。下面我们讨论实验引例中的单摆问题.练习3求方程0)0(',)0(,sin"0mgml的数值解.不妨取15)0(,8.9,1gl.则上面方程可化为0)0(',15)0(,sin8.9"先看看有没有解析解.运行MATLAB代码clear;s=dsolve('D2y=9.8*sin(y)','y(0)=15','Dy(0)=0','t')simplify(s)知原方程没有解析解.下面求数值解.令',21yy可将原方程化为如下方程组0)0(,15)0()sin(8.9''211221yyyyyy建立M文件fun9.m如下%M文件fun9.mfunctionf=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'),ylabel('y1')结果见图6.7.20123456789101515.51616.5ty1图6.7.2数值解由图6.7.2可见，随时间t周期变化。数学实验-常微分方程4/9练习4(刚性方程组求解)求下面刚性微分方程的解1)0(,2)0(,100'99.9901.0'2122211yyyyyyy使用dsolve可知解析解为)100exp(),100exp()01.0exp(21tytty下面求数值解.建立M文件fun10.m如下%M文件fun10.mfunctio...