MATLAB求解方程解析解和数值解

辽宁工程技术大学上机实验报告成绩实验名称Matlab求解方程的数值解和解析解院系专业班级姓名学号日期实验目的简述本次实验目的：1、熟悉MATLAB软件环境；2、熟悉MATLAB的常用运算符；3、了解MATLAB的一些常用函数；实验准备你为本次实验做了哪些准备：提前熟悉线性代数中的方程求解相关运算；提前熟悉Matlab中的方程求解相关的命令；实验进度本次共有4个练习，完成4个。实验总结日本次实验的收获、体会、经验、问题和教训：通过本次实验我发现，在Matlab中一些算法会变得很简单，有时候并不需要我们去了解具体的程序内部的算法，只要我们学会如何熟练运用Matlab软件就好。学会如何运用Matlab中的算法会对我们研究一些问题带来很大的方便，解决问题会变得很方便，免去了一些手动难以解决的问题。教师评语用MATLAB求解质点振动方程振动是日常生活和工程技术中常见的一种运动形式。利用常系数线性微分方程的理论来讨论有关自由振动和强迫振动的相关问题。利用MATLAB数学软件大致可分四类情况：（1）无阻尼自由振动情况；（2）有阻尼自由振动；（3）无阻尼强迫振动；（4）有阻尼强迫振动求其数值解和解析解；MATLAB软件求解微分方程解析解的命令“dsolve()”求通解的命令格式：（’微分方程’,’自变量’）注：微分方程在输入时，一阶导数y’应输入Dy，y’’应输入D2y等，D应大写。1，无阻尼自由振动情况：常见的数学摆的无阻尼微小振动方程代码如下：t=0:pi/50:2*pi;y=2*sin(3*t+2);plot(t,y,'b')01234567-2-1.5-1-0.500.511.522，有阻尼自由振动由无阻尼振动的通解可以看出，无阻尼振动是按照正弦规律运动的，摆动似乎可以无限期的进行下去，但事实上，空气从在阻力，在运动时，我们必须把空气阻力考虑在内，所以我们得到有阻尼摆动方程为：记u/m=2n,g/l=w^2,这里n,w是正常数，所以：y=dsolve('D2y+2*n*Dy+w^2*y=0','t');（4.43）解得:y=C3*exp(-t*(n+((n+w)*(n-w))^(1/2)))+C2*exp(-t*(n-((n+w)*(n-w))^(1/2)))（1）小阻尼情形:nw时，方程（4.43）的通解为：y=exp(-n*t)*(c1*cos(w1*t)+c2*sin(w1*t))和前面无阻尼的情形一样，可以把上式的通解改写为一下形式：y=A*exp(-n*t)*sin(w1*t+Q),(4.45)这里的A，Q为任意常数。用matlab操作得到：t=0:0.1:10;y=3*exp(-0.1*t).*sin(5*t+4);plot(t,y,'k-')如图：012345678910-3-2-10123由（4.45）可见，摆动的运动不是周期的，振动的幅度随着时间的增加而不断减小。（2）大阻尼情形：你w时；r2r10;方程（4.43）的通解：Y=c1*exp(r1*t)+c2*exp(r2*t),这里的c1,c2为任意常数；（3）临界的情形：即，n=w的情形，方程（4.43）的通解解为y=exp(-n*t)*(c1+c2*t),这里的，c1，c2为任意常数；由MAYLAB绘制图像得：t=0:0.1:100;yy=exp(-0.2*t).*(0.5-0.2*t);plot(t,yy,'r-')0102030405060708090100-0.3-0.2-0.100.10.20.30.40.53.有阻尼自由振动无阻尼自由振动和有阻尼自由振动都属于自由振动，它对应于一个二阶常数齐次线性微分方程。当一个振动系统还经常受到一个外力作用时，这种振动称为强迫振动。＝Ａｓｉｎ（ｔ＋）＋Ｈ／（２－ｐ２）ｓｉｎ（ｐｔ）取Ａ＝２；＝５；ｐ＝３；t=0:pi/50:2*pi;y=2*sin(5*t+2)+1/9*sin(3*t);plot(t,y,'k')01234567-2.5-2-1.5-1-0.500.511.522.54.有阻尼强迫振动这时摆动的运动方程（1.11）变为：D2y+2*n*Dy+w.^2*y=H*sin(p*t),（4.52）根据实际情况，我们只讨论小阻尼的情况，即nw的情形，这时的（4.52）对应的齐次线性方程的通解为：Y=A*exp(-n*t)*sin(w1*t+Q),(4.45)这里的A，Q为任意常数,w1=sqrt(w^2-n^2)现在求（4.52）的一个特解，这时可以寻求如：Y=M*cos(p*t)+N*sin(p*t)的特解，这里M，N是待定系数，代入（4,53），（4,52）化简得：Y=H*sin(Q)*cos(p*t)+H*cos(Q)*sin(P*t)=H*sin(p*t+Q),因此，（4.52）的通解为：Y=A*exp(-n*t)*sin(w1*t+Q)+H/sqrt((w^2-p^2)^2+4*n^2*p^2)*sin(p*t+Q);Matlab操作如下：A=3;n=0.2;w1=2;Q=0.5;H=5;w=3;p=2;t=0:0.1:50;Y=A*exp(-n*t).*sin(w1*t+Q)+H/sqrt((w.^2-p^2)^2+4.*n.^2.*p^2).*sin(p.*t+Q);plot(t,Y,'k-')05101520253035404550-3-2-101234利用MATLAB求解方程的数值解和解析解用龙格-库塔法求方程数值解：方程：i(t)=si-0.3ii(0)=0.02s(t)=-sis(0)=0.98建立函数M-文件：functiony=ill(t,x)a=1;b=0.3;y=[a*x(1)*x(2)-b*x(1),-a*x(1)*x(2)]';建立程序命令：ts=0:50;x0=[0.02,0.98];[t,x]=ode45('ill',ts,x0);[t,x]plot(t,x(:,1),t,x(:,2));grid,pause0510152025303540455000.10.20.30.40.50.60.70.80.91plot(x(:,2),x(:,1));grid00.10.20.30.40.50.60.70.80.9100.050.10.150.20.250.30.35用dsolve求方程的解析解：方程：x’=1-x^2x(0)=0clearx=dsolve('Dx=1-x^2','x(0)=0')x=tanh(t)方程：x’=1-x^2x(-3)=0.99clearx=dsolve('Dx=1-x^2','x(-3)=0.99')x=tanh(t+atanh(99/100)+3)