常微分方程实验报告

1、常微分方程课程实验报告实验名称 Matlab在常微分方程中的应用 班级信息1101学号201110010122姓名杨晓冰实验地点数学实验中心实验时间2013.6.7下午任课教师 孔令才评分一、 实验目的1. 掌握用Matlab求常微分方程（组）解析解的方法；2. 了解用Matlab求常微分方程（组）数值解的常用方法；3. 掌握Matlab作图方法；4. 培养编程与上机调试能力。二、用到的命令或函数 1. 常用作图函数 1.1 二维图形的绘制 plot, ezplot, 1.2 三维图形的绘制 plot3, ezplot3,mesh, meshgrid, surf;2. 求常微分方程（组）的解析

2、解函数 dsolve, 3. 求常微分方程（组）的数值解函数 ode23, ode45,4. 化简 simplify 三、实验内容 1. 求下列方程（组）的通解，并作出解的图形： （1） ； （2） ； （3） 2. 求下列方程（组）的特解，并作出解的图形： （1） ； （2） ； （3） 3. 求初值问题的数值解，求解范围为0,1，并做出图形。 4. 求初值问题的特解（精确解和近似解），求解范围为0,1，比较两种解的误差并作图观察。 5. 求微分方程组满足初始条件的特解（精确解），并画出解函数图形；再分别用ode23,ode45求此问题的数值解（近似解），求解区间为0,2, 画图比较两种解的

3、误差。实验步骤第1题（1）：【1】编写脚本文件chang1.msyms x yy=dsolve('Dy=y/x*(1+log(y)-log(x)','x')【2】在command窗口运行chang1.m,得到的结果：>> chang1 y = x/exp(C1*x)【3】作图：for C=0:0.01:1%当C1取大于0时hold on;x1=0:0.1:10;y1=x1./exp(C.*x1);subplot(1,2,1),plot(x1,y1),legend('C1>0')endfor c=-1:0.01:1 %当C1取小于

4、0时 hold on; x2=0:0.1:10; y2=x2./exp(c.*x2); subplot(1,2,2),plot(x2,y2),legend('C1<0')end【4】图像为：第1题（2）：【1】编写脚本M文件chang2.msyms x yy=dsolve('D2y-2*Dy+5*y-exp(x)*sin(2*x)','x')【2】在command窗口运行chang2.m,得到的结果：>> chang2 y = -1/8*exp(x)*(-8*sin(2*x)*C2-8*cos(2*x)*C1-sin(2*x)+

5、2*cos(2*x)*x)【3】作图：for c=0:1:10%C1和C2都大于 0 for c1=0:1:10 hold on; x1=40:0.1:50;%x取值4050 y1= -1/8.*exp(x1).*(-8.*sin(2.*x1)*c1-8.*cos(2.*x1)*c-sin(2.*x1)+2.*cos(2.*x1).*x1); subplot(2,2,1), plot(x1,y1),legend('C1>0,C2>0') endendfor c2=0:1:10%C1大于0，C2小于0 for c3=-10:1:0 hold on; x2=40:0.1

6、:50;%x取值4050 y2= -1/8.*exp(x2).*(-8.*sin(2.*x2)*c3-8.*cos(2.*x2)*c2-sin(2.*x2)+2.*cos(2.*x2).*x2); subplot(2,2,2), plot(x2,y2),legend('C1>0,C2<0') endend for c4=-10:1:0%C1小于0，C2大于0 for c5=0:1:10 hold on; x3=40:0.1:50;%x取值 4050 y3= -1/8.*exp(x3).*(-8.*sin(2.*x3)*c5-8.*cos(2.*x3)*c4-sin(

7、2.*x3)+2.*cos(2.*x3).*x3); subplot(2,2,3), plot(x3,y3),legend('C1<0,C2>0') endend for c6=-10:1:0%C1小于0，C2小于0 for c7=-10:1:0 hold on; x4=40:0.1:50;%x取值4050 y4= -1/8.*exp(x4).*(-8.*sin(2.*x4)*c7-8.*cos(2.*x4)*c6-sin(2.*x4)+2.*cos(2.*x4).*x4); subplot(2,2,4), plot(x4,y4),legend('C1<

8、;0,C2<0') endend 【4】图像为：第1题(3):【1】编写脚本M文件chang3.msyms x y z tx y z=dsolve('Dx=2*x-3*y+3*z','Dy=4*x-5*y+3*z','Dz=4*x-4*y+2*z','t');x=simplify(x)y=simplify(y)z=simplify(z)【2】在command窗口运行chang3.m,得到的结果：>> chang3 x = C1*exp(-t)+C3*exp(2*t) y = C1*exp(-t)+C2*e

9、xp(-2*t)+C3*exp(2*t) z = C2*exp(-2*t)+C3*exp(2*t)【31】作图(先做出x关于t，y关于t,z关于t的函数图像)：hold on;for c1=0:0.1:1%C1、C2和C3都大于0 for c2=0:0.1:1 for c3=0:0.1:1 t=0:0.1:1; x=c1.*exp(2.*t)+c3.*exp(-t); y=c1.*exp(2.*t)+c2.*exp(-t)+exp(-2.*t)*c3; z=c2.*exp(2.*t)+exp(-2.*t)*c3; subplot(1,3,1),plot(t,x),legend('tx&

10、#39;) %x关于t的函数图像 subplot(1,3,2),plot(t,y),legend('ty') % y关于t的函数图像 subplot(1,3,3),plot(t,z),legend('tz') % z关于t的函数图像 end endend【3.2】作图（再做出x、y、z的三维图像）：hold on;X=;Y=;Z=;for c1=0:0.2:1%C1、C2和C3都大于0 for c2=0:0.2:1 for c3=0:0.2:1 t=0:0.1:1; x=c1.*exp(2.*t)+c3.*exp(-t); y=c1.*exp(2.*t)+c2.

11、*exp(-t)+exp(-2.*t)*c3; z=c2.*exp(2.*t)+exp(-2.*t)*c3; X=X,x; Y=Y,y; Z=Z,z; end endend plot3(X,Y,Z,'y-','LineWidth',1) %x、y、z的三维立体图像【4.1】图像为(【3.1】对应的)：【4.2】图像为（【3.2】对应的）：第2题（1）：【1】编写脚本M文件 chang4.msyms x yy=dsolve('Dy+2*y-4*sin(100*pi*x)','y(0)=0','x')【2】在comma

12、nd窗口运行chang4.m,得到的结果为：>> chang4 y = 100*exp(-2*x)*pi/(1+2500*pi2)-2*(50*cos(100*pi*x)*pi-sin(100*pi*x)/(1+2500*pi2)【3】作图：hold on;Y1=;for x1=1:0.1:10 y1=100.*exp(-2.*x1)*pi/(1+2500*pi2)-2*(50.*cos(100*pi.*x1)*pi-sin(100*pi.*x1)/(1+2500*pi2) Y1=Y1,y1;endx1=1:0.1:10plot(x1,Y1,'r-')【4】图像为：

13、第2题（2）：【1】编写脚本M文件chang5.msyms x y y=dsolve('D2y+4*Dy+29*y=0','y(0)=0,Dy(0)=15')【2】在command窗口运行chang5.m,得到的结果为：>> chang5 y = 3*exp(-2*t)*sin(5*t)【3】作图：hold on;Y1=;for t=0:0.01:10 y1=3.*exp(-2.*t).*sin(5.*t); Y1=Y1,y1;end t=0:0.01:10 plot(t,Y1,'y*')【4】图像为：第2题（3）：【1】编写脚本M文

14、件chang6.mx,y=dsolve('Dx=-3*x+y','Dy=8*x-y','x(0)=1','y(0)=1','t');x=simplify(x)y=simplify(y)【2】在command窗口运行chang6.m,得到的结果为：>> chang6 x = exp(t) y = 4*exp(t)【31】作图（先做出x关于t,y关于t的函数图像）：for t=0:0.1:10 x1=exp(t);subplot(1,2,1),plot(t,x1,'y*'),legend(&

15、#39;tx图像') hold on y1=4.*exp(t); subplot(1,2,2),plot(t,y1,'g*'),legend('ty图像')end【3.2】作图（再做出x关于t、y关于t合在一起的函数图像）： X1=; Y1=;for t=0:0.1:10 x1=exp(t); y1=4.*exp(t); hold on; X1=X1,x1; Y1=Y1,y1;endt=0:0.1:10;plot( t,X1,'y*',t,Y1,'g*')【3.3】作图（最后做出x、y的二维图像）： X1=; Y1=;fo

16、r t=0:0.1:10 x1=exp(t); y1=4.*exp(t); hold on; X1=X1,x1 Y1=Y1,y1endplot(X1,Y1,'b-','LineWidth',4)xlabel('x')ylabel('y')【4.1】图像为（【3.1】对应的）：【4.2】图像为（【3.2】对应的）：【4.3】图像为（【3.3】对应的）：第3题：【1】编写函数文件cwf.mfunction dy=cwf(x,y)dy=-2*y+2*x2+2*x;【2】编写脚本M文件chang7.mX,Y=ode45('cwf&

17、#39;,0,1,1)plot(X,Y,'r-','LineWidth',2)【3】图像为：【4】数值解的取值为（当x在01区间内取值时，ode45的方法）：>> Y'= Columns 1 through 9 1.0000 0.9519 0.9073 0.8663 0.8287 0.7944 0.7633 0.7353 0.7103 Columns 10 through 18 0.6883 0.6690 0.6526 0.6388 0.6277 0.6191 0.6130 0.6093 0.6080 Columns 19 through 2

18、7 0.6091 0.6124 0.6179 0.6256 0.6354 0.6473 0.6612 0.6771 0.6950 Columns 28 through 36 0.7149 0.7366 0.7602 0.7856 0.8129 0.8419 0.8727 0.9052 0.9394 Columns 37 through 410.9753 1.0129 1.0521 1.0929 1.1353第4题：【1】编写函数M文件cwf1.m(用于求近似解)function dy=cwf1(x,y)dy=(cos(x)-2*x*y)/(x2-1)【2】编写脚本M文件chang8.mX,Y=

19、ode45('cwf1',0,1,1)y1=Y %转置【3】则求出的近似解为（x取值为01）：>>y1= Columns 1 through 10 1.0000 0.9756 0.9524 0.9303 0.9093 0.8892 0.8701 0.8520 0.8347 0.8183 Columns 11 through 20 0.8028 0.7880 0.7742 0.7611 0.7488 0.7374 0.7269 0.7172 0.7085 0.7008 Columns 21 through 30 0.6941 0.6886 0.6843 0.6815

20、0.6802 0.6809 0.6837 0.6890 0.6976 0.7102 Columns 31 through 40 0.7277 0.7519 0.7851 0.8319 0.8964 0.9910 1.1404 NaN NaN NaN Column 41 -Inf【4】在command窗口运行以下语句，用于求精确解：>> y=dsolve('(x2-1)*Dy+2*x*y-cos(x)=0','y(0)=1','x') %精确解 y = (sin(x)-1)/(x2-1)【5】作图：X,Y=ode45('cwf1

21、',0,1,1)Y1=;for x=0:0.1:1 y1=(sin(x)-1)./(x2-1) Y1=Y1,y1;endx=0:0.1:1plot(x,Y1,'g-','LineWidth',2)hold onplot(X,Y,'y-','LineWidth',2),legend('jin si jie')第5题：(1)求精确解：【1】编写脚本M文件chang9.m(用于求精确解)x,y=dsolve('Dx+x+y=0','Dy+x-y=0','x(0)=1'

22、;,'y(0)=0','t')【2】在command窗口运行chang9.m，得到的结果为（精确解）：>> chang9 x = (-1/4*2(1/2)+1/2)*exp(2(1/2)*t)+(1/4*2(1/2)+1/2)*exp(-2(1/2)*t) y = -(-1/4*2(1/2)+1/2)*2(1/2)*exp(2(1/2)*t)+(1/4*2(1/2)+1/2)*2(1/2)*exp(-2(1/2)*t)-(-1/4*2(1/2)+1/2)*exp(2(1/2)*t)-(1/4*2(1/2)+1/2)*exp(-2(1/2)*t)【3】

23、作图：X1=;Y1=;for t=0:0.1:2 hold on; x1=(-1/4*2(1/2)+1/2)*exp(2(1/2)*t)+(1/4*2(1/2)+1/2)*exp(-2(1/2)*t); y1=-(-1/4*2(1/2)+1/2)*2(1/2)*exp(2(1/2)*t)+(1/4*2(1/2)+1/2)*2(1/2)*exp(-2(1/2)*t)-(-1/4*2(1/2)+1/2)*exp(2(1/2)*t)-(1/4*2(1/2)+1/2)*exp(-2(1/2)*t); X1=X1,x1; Y1=Y1,y1;endt=0:0.1:2 plot(t,X1,'y-&#

24、39;,'LineWidth',2) plot(t,Y1,'r-','LineWidth',2)【4】图像为：（2）求近似解（用ode45求近似解）：【1】编写函数M文件cwf2.mfunction dy=cwf2(t,y)dy=zeros(2,1);dy(1)=-y(1)-y(2);dy(2)=-y(1)+y(2)【2】编写脚本M文件chang10.mt,y=ode45('cwf2',0,2,1,0);plot(t,y(:,1),'b d',t,y(:,2),'g*')y1=y'%转置【3

25、】图像为：【4】用ode45求得的近似解为（x取值为02）：y1 = Columns 1 through 10 1.0000 0.9999 0.9999 0.9998 0.9998 0.9995 0.9993 0.9990 0.9988 0.9975 0 -0.0001 -0.0001 -0.0002 -0.0002 -0.0005 -0.0007 -0.0010 -0.0012 -0.0025 Columns 11 through 20 0.9963 0.9951 0.9938 0.9876 0.9816 0.9756 0.9696 0.9411 0.9145 0.8896 -0.0037

26、-0.0050 -0.0062 -0.0125 -0.0188 -0.0251 -0.0314 -0.0628 -0.0944 -0.1262 Columns 21 through 30 0.8665 0.8332 0.8041 0.7790 0.7578 0.7404 0.7267 0.7167 0.7102 0.7073 -0.1582 -0.2099 -0.2626 -0.3167 -0.3723 -0.4298 -0.4894 -0.5515 -0.6163 -0.6843 Columns 31 through 40 0.7079 0.7121 0.7198 0.7311 0.7461

27、 0.7648 0.7873 0.8138 0.8443 0.8791 -0.7556 -0.8307 -0.9100 -0.9939 -1.0827 -1.1769 -1.2770 -1.3835 -1.4969 -1.6178 Columns 41 through 50 0.9183 0.9620 1.0106 1.0642 1.1232 1.1877 1.2583 1.3351 1.4185 1.5091 -1.7468 -1.8845 -2.0317 -2.1890 -2.3573 -2.5374 -2.7301 -2.9365 -3.1576 -3.3945 Columns 51 t

28、hrough 60 1.6072 1.7134 1.8281 1.9520 2.0857 2.2297 2.3850 2.4199 2.4554 2.4915 -3.6484 -3.9205 -4.2123 -4.5251 -4.8605 -5.2203 -5.6061 -5.6928 -5.7808 -5.8701 Column 61 2.5282 -5.9608(3)求近似解（用ode23求近似解）：【1】编写函数M文件cwf2.mfunction dy=cwf2(t,y)dy=zeros(2,1);dy(1)=-y(1)-y(2);dy(2)=-y(1)+y(2)【2】编写脚本M文件ch

29、ang10.mt,y=ode23('cwf2',0,2,1,0);plot(t,y(:,1),'b d',t,y(:,2),'g*')y1=y'%转置【3】图像为：【4】用ode23求得的近似解为（x取值为02）：y1 = Columns 1 through 10 1.0000 0.9999 0.9995 0.9975 0.9877 0.9413 0.8685 0.7933 0.7346 0.7075 0 -0.0001 -0.0005 -0.0025 -0.0125 -0.0626 -0.1553 -0.2846 -0.4525 -0.6672 Colu