常微分方程實(shí)驗(yàn)報(bào)告

1、常微分方程課程實(shí)驗(yàn)報(bào)告實(shí)驗(yàn)名稱 Matlab在常微分方程中的應(yīng)用 班級(jí)信息1101學(xué)號(hào)201110010122姓名楊曉冰實(shí)驗(yàn)地點(diǎn)數(shù)學(xué)實(shí)驗(yàn)中心實(shí)驗(yàn)時(shí)間2013.6.7下午任課教師 孔令才評(píng)分一、 實(shí)驗(yàn)?zāi)康?. 掌握用Matlab求常微分方程（組）解析解的方法；2. 了解用Matlab求常微分方程（組）數(shù)值解的常用方法；3. 掌握Matlab作圖方法；4. 培養(yǎng)編程與上機(jī)調(diào)試能力。二、用到的命令或函數(shù) 1. 常用作圖函數(shù) 1.1 二維圖形的繪制 plot, ezplot, 1.2 三維圖形的繪制 plot3, ezplot3,mesh, meshgrid, surf;2. 求常微分方程（組）的解析

2、解函數(shù) dsolve, 3. 求常微分方程（組）的數(shù)值解函數(shù) ode23, ode45,4. 化簡(jiǎn) simplify 三、實(shí)驗(yàn)內(nèi)容 1. 求下列方程（組）的通解，并作出解的圖形： （1） ； （2） ； （3） 2. 求下列方程（組）的特解，并作出解的圖形： （1） ； （2） ； （3） 3. 求初值問題的數(shù)值解，求解范圍為0,1，并做出圖形。 4. 求初值問題的特解（精確解和近似解），求解范圍為0,1，比較兩種解的誤差并作圖觀察。 5. 求微分方程組滿足初始條件的特解（精確解），并畫出解函數(shù)圖形；再分別用ode23,ode45求此問題的數(shù)值解（近似解），求解區(qū)間為0,2, 畫圖比較兩種解的

3、誤差。實(shí)驗(yàn)步驟第1題（1）：【1】編寫腳本文件chang1.msyms x yy=dsolve('Dy=y/x*(1+log(y)-log(x)','x')【2】在command窗口運(yùn)行chang1.m,得到的結(jié)果：>> chang1 y = x/exp(C1*x)【3】作圖：for C=0:0.01:1%當(dāng)C1取大于0時(shí)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 %當(dāng)C1取小于

4、0時(shí) 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窗口運(yùn)行chang2.m,得到的結(jié)果：>> 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窗口運(yùn)行chang3.m,得到的結(jié)果：>> 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關(guān)于t，y關(guān)于t,z關(guān)于t的函數(shù)圖像)：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關(guān)于t的函數(shù)圖像 subplot(1,3,2),plot(t,y),legend('ty') % y關(guān)于t的函數(shù)圖像 subplot(1,3,3),plot(t,z),legend('tz') % z關(guān)于t的函數(shù)圖像 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】對(duì)應(yīng)的)：【4.2】圖像為（【3.2】對(duì)應(yīng)的）：第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窗口運(yùn)行chang4.m,得到的結(jié)果為：>> 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窗口運(yùn)行chang5.m,得到的結(jié)果為：>> 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窗口運(yùn)行chang6.m,得到的結(jié)果為：>> chang6 x = exp(t) y = 4*exp(t)【31】作圖（先做出x關(guān)于t,y關(guān)于t的函數(shù)圖像）：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關(guān)于t、y關(guān)于t合在一起的函數(shù)圖像）： 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】對(duì)應(yīng)的）：【4.2】圖像為（【3.2】對(duì)應(yīng)的）：【4.3】圖像為（【3.3】對(duì)應(yīng)的）：第3題：【1】編寫函數(shù)文件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】數(shù)值解的取值為（當(dāng)x在01區(qū)間內(nèi)取值時(shí)，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】編寫函數(shù)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 %轉(zhuǎn)置【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ùn)行以下語(yǔ)句，用于求精確解：>> 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窗口運(yùn)行chang9.m，得到的結(jié)果為（精確解）：>> 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】編寫函數(shù)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'%轉(zhuǎn)置【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】編寫函數(shù)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'%轉(zhuǎn)置【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