MATLAB实验指导书加程序+上机实例

1、MATLA语言实验指导书华东交通大学电气学院张永贤2006年2月实验一 MATLAB工作环境熟悉及简单命令的执行一、 实验目的：熟悉MATLAB的工作环境，学会使用MATLAB进行一些简单的运算。二、实验内容：MATLAB勺启动和退出，熟悉MATLAB勺桌面(Desktop )，包括菜单(Menu)、 工具条 (Toolbar )、命令窗口 (CommancWindow)、历史命令窗口、工作空间(Workspace) 等；完成一些基本的矩阵操作；学习使用在线帮助系统。三、实验步骤：1 启动 MATLAB熟悉 MATLA啲桌面。2、在命令窗口执行命令完成以下运算，观察workspace的变化，记

2、录运算结果。(1) (365-52 2-70 )3(2) >>area=pi*2.5A2(3、已知 x=3, y=4,在 MATLAB中求 z :23z xy2x y(4、将下面的矩阵赋值给变量ml,在workspace中察看ml在内存中占用的字节数。162313511108m1 =97612414151执行以下命令>>m1( 2,3 )>>m1( 11 )>>m1( : , 3 )>>m1( 2 : 3,1 : 3 )>>m1( 1 ,4 ) + m1( 2 ,3 ) + m1( 3 ,2 ) + m1( 4 ,1)(5

3、、执行命令>>help abs查看函数abs的用法及用途，计算 abs( 3 + 4i )(6、执行命令>>x=0:0.1:6*pi;>>y=5*s in (x);>>plot(x,y)(6、运行MATLAB勺演示程序，>>demo以便对 MATLAB有一个总体了解。四、思考题1、以下变量名是否合法？为什么？(1) x2(2) 3col(3) _row(4) for2、求以下变量的值，并在 MATLAB验证。(1) a = 1 : 2 : 5 ;(2) b = a' a' a'(3) c = a + b ( 2

4、 ,:)实验二MATLAB语言矩阵运算、实验目的：掌握基本的矩阵运算及常用的函数。、实验内容：1 232141abc04 5613 521、下列运算是否合法，为什么？如合法，结果是多少? resultl = a'(2) result2 = a * b(3) result3 = a + b(4) result4 = b * d(5) result5 = b ; c' * d(6) result6 = a . * b(7) result7 = a . / b(8) result8 = a . * c(9) result9 = a . b(10) result10 = a . A2(

5、11) resultll = a 人2(12) resultll = 2 . a a2、用MATLAB求下面的的方程组。7212X1491532X27(1)22115X3113213X40xy z1x2yz w8(2)2xy3w33x3y5z6w57212915323、已知A2211513213(1)求矩阵A的;秩(ranik)(2)求矩阵A的;行列式:(detenminant)(3)求矩阵A的:逆 (inverse)(4)求矩阵A的!持征值【及特征向量(eigien value and eige nvector)4、关系运算与逻辑运算已知 a=20,b=-2,c=0,d=1(1) r1 =

6、a > b(2) r2 = a > b & c > d(3) r3 = a = b* (-10) r4 = b | c、思考题102门21029n 1091022,求 y= ?(用 format long 查看 y 的值)实验三程序的编辑及调试MATLAB序编辑、运行及调试方法。、实验目的：掌握 、实验内容：启动MATLAB后，(Editor/Debugger点击File|New|M-File ，启动 MATLAB的程序 编辑及 调试器 )，编辑以下程序，点击File|Save 保存程序，注意文件名最好用英 文字符。点击 Debug|Run运行程序，在命令窗口查看运行

7、结果，程序如有错误则改正。注：数论中一个有趣的题目：任意一个正整数，若为偶数，则用2除之，若为奇数，则与3相乘再加上236运行下面的程序，按程序提示输入n=123,5,7等数来验证这一结论。%classic "3n+1" problem from n umber theory.while 1若为偶数，则用重复此过程，最终得到的结果为 1。如：110 53 1016168 4 21n=in put('En ter n,n egative quits:'if n<=0)；breakenda=n;while n>1if rem(n,2)=0 n=n /

8、2;elsen=3*n+1;enda=a ,n ;endaend2、编程求满足m ii12i10000的最小m直。三、思考题用对分法求解方程2e xsinx在0, 1内的解，并验证，在程序中统计出对分次数。提示：先将原方程转化成f(x) 2ex sin x0的形式。对分法的基本思想是：一个一元方程有实数解。取该区间的中点xm=(x1+x2)/2f(x)=0，若，判定f(x1)f(x1)*f(x2)<0 ，则在x1,x2区间内 和f(x2)二者中哪一个与f(xm)异号,若f(x1)*f(xm)<0 ，则解存在的区间缩小为x1,xm，否则解存在的区间缩小为xm,x2。重复这样的步骤，直

9、到区间的长度小于一个可以接受的小数(比如1e-10 ),则认为中点即是原方程的解。实验四 函数的编写及调试一、实验目的：掌握MATLAB函数的编写及调试方法。二、实验内容：1、编写一个函数，计算下面函数的值，给出标量x的值，调用该函数后，返回y的值。function y=myfun1(x)sin x,x 0y(x) x,0x3x 6,x 3选择一些数据测试你编写的函数。2、编写一个函数求向量 x中元素的平均值、最大值、最小值、均方根值。function m_x,max_x,min_x,rms_x=myfun2(x)方均根值(RootMean Square)的计算公式为：rms用下面数据测试你写

10、的函数：(1) x=sin(0:0.01:6*pi)(2) x=rand(1,200)，得到的x为200个(0, 1)之间均匀分布的随机数。3、编写一个函数，给出一个向量xxpX2,xn，生成如下范德蒙矩阵。function v=myvander(x)111X1X2Xn222X1X2Xnn 1n 1n 1X1X2xn例如：>>v=myvander(2 3 4 5)得v=11112345491625827 64 125生成一些数据测试你写的函数。三、思考题编写程序，用如下迭代公式求 a，a的值分别为：3,17，113。迭代的终止条件为Xn1 Xn10，迭代初值X。1.0，迭代次数不超

11、过100次。分别对迭代结果和迭 代 次 数。x2a 2x2x2XnXn 122xn 1i、实验五在同一坐标系下绘制下面三个函数在tMATLAB的绘图0，4 的图象。2、y 3编写程序，t、t4 e 0.1t sin( t )绘制下面函数在区间-6 , 6中的图象。选择合适的步距,sin x,y(x)x,x 6,3、用compass函数画下面相量图ua =1; ub =uc=cos(2*pi/3)+si n(2*pi/3)*i;compass(ua,ub,uc,ua-ub,ub-uc,uc-ua)cos(-2*pi/3)+si n(-2*pi/3)*i4、三维空间曲线绘制 z=0:0.1:4*p

12、i; x=cos(z);y=s nz）；plot3（x,y,z）x和y的取值范围设为-3 ，5、用mesh或surf函数，绘制下面方程所表示的三维空间曲面,3。10 10三、思考题t -2，2 范围内的图象。在同一坐标系下，用不同颜色和线型绘制以下两个函数在y12°刖y22e0.21实验六MATLAB数值运算一、实验目的：掌握MATLAB常用的数值运算函数。二、实验内容：54321、求代数方程3x 4x 7x 2x 9x 120的5个根，并将其用星号（*）标记在复平面图上。（用roots和plot函数）。5 丿2、 求代数方程x 1 0的5个根，并将其用星号（* ）标记在复平面图上。

13、（用roots和 plot函数）。3、求下面函数在0.5,4区间内的过零点。（用fzero函）f(x) x322x sin(x) 5xcos(x)4、RQ»i +50DU3Lo已知R=50欧姆，U=4V二极管D正向电流与电压的关系为:u dqKTIse其中：Ud为二极管正向电压Is为反向饱合电流，取 10-12AK为玻尔茨曼常数，1.38*10 -23T为绝对温度，取 300开尔文（27摄氏度）q为电子电荷1.6*10 C求此电路中的电流Id和二极管正向电压 Ud (要求用fsolve函数求解)5、实验数据处理：已知某压力传感器的测试数据如下表p0.01.12.12.84.25.06

14、.16.98.19.09.9u101113141718222429343932p为压力值，u为电压值，试用多项式u(p) ap bp cp d来拟合其特性函数，求出a,b,c,d ，并把拟合曲线和各个测试数据点画在同一幅图上。实验七MATLAB用1、以原点为奇对称中心的方波y(wt)，可以用相应频率的基波及其奇次谐波合成。411y(wt)sin wtsin3wtsin 5wt35n 1,2,3,1(2n 1)sin(2n1)wt取的阶数越多，越接近方波，但总消除不了边缘上的尖峰，这称为吉布斯效应。设方波频率为50Hz,时间 t 取00.04 秒(f=50Hz,w=2*pi*f,h=1e-5,t

15、f=40e-3,t= 0:h:tf)，编写程序，画出如下用1次谐波、1,3次谐波、1,3,5,7,9次谐波，1,3,5,19次谐波合成的近似方波。(产 生方波的函数为：square)2、用Simulink求解下图所示电路0100微秒内的响应。已知R=6*104欧，C=1700 微法，L=6*10-9 享，uc(0)=15kV。Gain2tr ScopeUcUcScope1模块参数设置：Integrator1的 Initial condition:15kV在命令窗口为R,L,C赋值。仿真参数设置如下：Start time:0Stop time:100e-6Solver Type:Variable

16、-step Solver:ode45Max step size:1e-7Min step size:autoIn itial step size:autoRelative tolera nce:1e-3Absolute tolera nce:1e-6MATLAB实验程序 实验1第1题x=2*si n( 85*pi/180”(1+exp(2)x =0.2375.x=2 1+2i;-0.45 5;y=0.5*log(x+sqrt(1+xA (2)y =0.7114 - 0.0253i 0.8968 + 0.3658i0.2139 + 0.9343i 1.1541 - 0.0044i 或x=2 1+

17、2i;-0.45 5;d=0.5*log(x+sqrt(1+x*x)d =0.7114 - 0.0253i 0.8968 + 0.3658i0.2139 + 0.9343i 1.1541 - 0.0044i 或x=2 1+2*i;-0.45 5;y=0.5*log(x+sqrt(1+xA (2)y =0.7114 - 0.0253i 0.8968 + 0.3658i0.2139 + 0.9343i 1.1541 - 0.0044i .a=-3.0:0.1:3.0;g=(exp(0.3*a)-exp(-0.3*a).*si n(a+0.3)/2+log(0.3+a)/2) 结果略>>

18、 t=0:0.5:2.5;>> f1=t.A2;>> f2=t.A2-1;>> f3=t.A2-2*t+1;>> z=(t>=0&t<1).*f1+(t>=1 &t<2).*f2+(t>=2 &t<3).*f3 z =00.250001.25001.0000 2.2500第2题(1)>> A=12 34 -4;34 7 87;3 65 7;>> B=1 3 -1;2 0 3;3 -2 7;>> A+6*Bans =1852-10467105215349

19、>> A-B+eye (3)ans =1232-2338851681(2)>> A*Bans =684462309-72596154-5241>> A.*Bans =1210246802619-13049>> AA3ans =372262338244860424737014918860076678688454142118820>> A.A3ans =172839304-643930434365850327274625343>> A/Bans =16.4000 -13.6000 7.600035.8000 -76.2000

20、50.200067.0000 -134.0000 68.0000 >> BAans =109.4000 -131.2000 322.8000 -53.0000 85.0000 -171.0000 -61.6000 89.8000 -186.2000>> A,Bans =1234-413-13478720336573-27>> A(1,3,:);BA2ans =1234-436574511101920-540第3题>> A=1 2 3 4 5;6 7 8 9 10;11 12 13 14 15;16 17 18 19 20;21 22 23 24

21、25;>> B=3 0 16;17 -6 9;0 23 -4;9 7 0;4 13 11;>> C=A*BC =9315077258335237423520397588705557753890717>> D=C(3:5,2:3)D =520397705557890717第4题(1)>> a=100:999;>> k=find(rem(a,21)=0);%找出能杯21整除的元素位置,find()函数找出不为 0的元素位置>> x=length(k)%获得向量k的元素个数并赋值给变量xx =43>> k=find(

22、rem(a,21)=0) % 显示能杯21整除的元素位置k =Colu mns 1 through 24627 486990 111 132 153 174 195 216 237 258279 300 321 342 363 384 405 426 447 468 489Columns 25 through 43510 531 552 573 594 615 636 657 678 699 720 741 762783 804 825 846 867 888>> y=100+k-1%显示能杯21整除的元素y =Colu mns 1 through 23105 126 147 16

23、8 189 210 231 252 273 294 315 336 357378 399 420 441 462 483 504 525 546 567Columns 24 through 43588 609 630 651 672 693 714 735 756 777 798 819 840861 882 903 924 945 966 987sh='CDe345Efg69K：>> k=fi nd(sh>='A'&sh<='Z');%找出大写字母的位置>> sh(k)=;%删除大写字母>> x=

24、sh(1:e nd)%显示处理后的字符x =e345fg69实验2第1题a=1 2 3;4 5 6;b=2 4 -1;1 3 5;c=1;0;-2;% a的转置d=1 4 7;8 5 2;3 6 0; >> result1=a' result1 =1 42 53 6>> result2=a*b >> result3=a+b result3 =362%error应采用点乘%求两个矩阵的和%巨阵相乘5 8 11 >> result4=b*d result4 =312222404913>> result5=b;c'*dres

25、ult5=312222404913-5-87>> result6=a.*bresult6 =28 -341530>> result7=a./bresult7 =0.50000.5000 -3.00004.00001.66671.2000>> result8=a.*c >> result9=a.b result9 =2.00002.0000 -0.33330.2500 0.60000.8333>> result10=a.A2 result10=1 4916 25 36 >> result1 仁aA2 >> res

26、ult12=2.Aa result12 =2 48163264第2题%矩阵点乘%巨阵右点除%error%矩阵左点除%errora和c维数不同等价于a*a(1)>> A=7 2 1 -2;9 15 3 -2;-2 -2 11 5;1 3 2 13;A/B 等价于 A*inv(B)>> B=4;7;-1;0;>> X=inv(A)*B%AB等价于 inv(A)*B,X =0.49790.14450.0629-0.0813 >> X1=ABX1 =0.49790.14450.0629-0.0813>> a=1 1 1 0;1 2 1 -1;

27、2 -1 0 -3;3 3 5 -6;>> b=1;8;3;5;>> x=in v(a)*bx =1.00005.0000-5.0000-2.0000第3题>> A=7 2 1 -2;9 15 3 -2;-2 -2 11 5;1 3 2 13;>> a1=ra nk(A)a1 =4>> a2=det(A)a2 =12568>> a3=i nv(A)a3 =0.1744-0.0303-0.01250.0270-0.10500.0789-0.01210.00060.00830.01730.0911-0.03110.0095-0

28、.0185-0.01030.0795>> V,D=eig(A)%为向量A的特征向量，D为特征值V =-0.76290.0919 + 0.0640i 0.0919 - 0.0640i -0.02990.62230.6087 + 0.0276i 0.6087 - 0.0276i 0.26370.0807-0.7474-0.74740.6434-0.15540.0342 - 0.2374i 0.0342 + 0.2374i 0.7180D =4.8554000012.6460 + 1.8333i000012.6460 - 1.8333i000015.8526第4题>> a=2

29、0;>> b=-2;>> c=0;>> d=1;>> r1=a>br1 =1>> r2=a>b&c>dr2 =0>> r3 = a = b* (-10)r3 =1>> r4 = b | cr4 =05>> y=0;for k=-10:10y=y+pow2(k);endformat lo ng>> yy =2.047999023437500e+003实验3方法一a=0;for i=1:20a=a+pow2(i);if a>10000m=i;breakende

30、ndm方法二a=0;i=1;while (a<10000)a=a+pow2(i);i=i+1;endm=i-1;m实验4第1题function y=myfun1(x)if x<=0y=si n( x);elseif x>0 &x<=3y=x;elseif x>3y=-x+6;end运行结果：>> y=myfu n1(-pi/2)y =-1>> y=myfu n1(0)y =0>> y=myfu n1(2)y =2>> y=myfun 1(4)y =2第2题fun ctio n m_x,max_x,min_x,

31、rms_x=myfu n2(x) %求平均值_sum_x=sum(x);%向量元素求禾口m,n=size(x);%最好用 n=length(x);m_x=sum_x/ n;%求最大值 采用逐个比较方式if x(1)>x(2)max_x=x(1);elsemax_x=x(2);endfor k=3:nif max_x<x(k)max_x=x(k);elsemax_x=max_x;%可省略endend%求最小值if x(1)<x(2)min_x=x(1);elsemin_x=x(2);endfor k=3:nif min_x>x(k) min_x=x(k);elsemin_

32、x=min_x;%可省略endend%求均方根值sum_x2=0;for k=1: nsum_x2=sum_x2+x(k).A2;rms_x=sqrt(sum_x2/n);endm_x;max_x;min_x;rms_x;%按照函数值行参顺序输出结果运行结果：>> m_x,max_x,min_x,rms_x=m yfun 2(si n(0:0.01:6*pi) m_x =-1.1256e-007max_x =1.0000min_x =-10000rms_x =0071>> m_x,max_x,min_x,rms_x=myfu n2(ra nd(1,200) m_x =&

33、quot;0.4977max_x =0.9961min_x =0046rms_x =0.5778第3题function v=myva nder(x)v1=va nder(x);%生成范德蒙矩阵v2=v1'v=flipud(v2);%实现矩阵上下翻转运行结果：>> v=myva nder(2:5)v =1111234549162582764125思考题function x,n=sqrt_a(a) x=1.0;for k=1:100m= x;x=x/2+a/(2*x);if abs(x-m)<=10A(-5) breakendendx;n=k;s=(x-sqrt(a);i

34、f s<=10A(-5)disp('正确');elsedisp('错误');end运行结果：>> x,n=sqrt_a(3)正确-x =1.7321n =5>> x, n=sqrt_a(17)正确-x =4.1231n =6>> x, n=sqrt_a(113)正确-x =10.6301n =8实验5第1题.方法1>> t=li nspace(0,4*pi,200); y1=t;y2=sqrt(t);y3=4*pi*exp(-0.1*t).*si n(t); plot(t,y1,'b',t,y

35、2,'g',t,y3,'r')方法2>> t=li nspace(0,4*pi,200); y1=t;y2=sqrt(t);y3=4*pi*exp(-0.1*t).*si n(t); t=t,t,t;y=y1,y2,y3;plot(t,y)第2题>> x=li nspace(-6,6,100); y=;for x0=xif x0<=0y=y,s in (x0);elseif x0<=3y=y,x0； elsey=y,-x0+6; endendplot(x,y);axis(-7 7 -2 4); title(' 分段函数

36、曲线'); text(-3*pi/2,1,'y=si n(x)'); text(2,2,'y=x');text(4,2,'y=-x+6');第3题ua=1;ub=cos(-2*pi/3)+si n(-2*pi/3)*i; uc=cos(2*pi/3)+si n(2*pi/3)*i; compass(ua,ub,uc,ua-ub,ub-uc,uc-ua) title('相量图);第4题>> z=0:0.1:4*pi; x=cos(z); y=s in( z); plot3(x,y,z,'rp'); tit

37、le('三维空间曲线');text(0,0,0,'origi n');xlabel( X),ylabel('Y'),zlabel('Z');grid;第5题(1)>>x=-3:0.1:3;x,y=meshgrid(x); z=-x.A2/10+y.A2/10; mesh(x,y,z);xlabel( X),ylabel('Y'),zlabel('Z'); title(' 立体网状图');(2)>>x=-3:0.1:3;x,y=meshgrid(x); z=-x

38、.A2/10+y.A2/10; surf(x,y,z);xlabel( X),ylabel('Y'),zlabel('Z'); title(' 立体曲面图');思考题>> t=-2*pi:pi/100:2*pi; y1=2.A(0.5*abs(t); y2=2*exp(-0.2*t); plot(t,y1,'-g');hold on;plot(t,y2,':r');legend('曲线 y1','曲线 y2'); hold off;grid;实验6第1题>>

39、A=3,4,7,2,9,12;x=roots(A)plot(x,'*');grid;x =-0.8612 + 1.4377i-0.8612 - 1.4377i0.6737 + 1.0159i0.6737 - 1.0159i-0.9583第2题>> A=1,0,0,0,0,-1;x=roots(A)plot(x,'*');grid;x =-0.8090 + 0.5878i-0.8090 - 0.5878i0.3090 + 0.9511i0.3090 - 0.9511i1.0000第3题%古计零点fplot('xA3+1/x',0.5,4

40、);hold on;fplot('2*xA2*si n(x)-5*x*cos(x)',0.5,4);hold off;m, n=gi nput%建立函数function y=f(x)y=xA3-2*xA2*si n(x)+5*x*cos(x)+1/x;%调用函数>> y仁 fzero('fz',1.5)y1 =1.5117>> y2=fzero('fz',2.5)y2 =2.6095第4题%古计零点axis(0,1,0,0.1);fplot('10A(-12)*exp(Ud*1.6*10A(-19)/(1.38*1

41、0A(-23)*300)-1) ',0,4); hold on;fplot('(Ud-4)/50',0,4);hold off;m, n=gi nput%建立函数fun ctio n f=myfu n(X)Id=X(1);Ud=X (2);f(1)=Id-10A(-12)*exp(Ud*1.6*10A(-19)/(1.38*10A(-23)*300)-1); f(2)=50*Id-Ud-4;%调用函数>> x=fsolve('myfu n',0,0.8,optimset('Display','off)x =0.0936

42、 0.6795%佥证结果>> K=m yfun(x)K =1.0e-008 *-0.30870第5题>> p=0.0,1.1,2.1,2.8,4.2,5.0,6.1,6.9,8.1,9.0,9.9; u=10,11,13,14,17,18,22,24,29,34,39;x=polyfit(p,u,3)%得多项式系数t=li nspace(0,10,100);y=polyval(x,t);%求多项式得值plot(p,u,'*',t,y,'r')%画拟和曲线x =0.0195 -0.04121.44699.8267实验7第1题>>

43、f=50;w=2*pi*f;h=1e-5;tf=40e-3;t=0:h:tf;wt=w*t;y1=4/pi*s in( wt);y2=4/pi*(si n(wt)+1/3*si n(3*wt);y3=4/pi*(si n(wt)+1/3*si n(3*wt)+1/5*si n(5*wt)+1/7*si n(7*wt)+1/9*si n(9*wt);y4=4/pi*(si n(wt)+1/3*si n(3*wt)+1/5*si n(5*wt)+1/7*si n(7*wt)+1/9*si n(9*wt)+1/11*s in( 11*wt)+1/13*si n(13*wt)+1/15*si n(15

44、*wt)+1/17*si n(17*wt)+1/19*si n(19*wt); y=square(wt);subplot(2,2,1);plot(wt,y1,wt,y);title('1 次谐波');subplot(2,2,2);plot(wt,y2,wt,y);title('1,3 次谐波');subplot(2,2,3);plot(wt,y3,wt,y);title('1,3,5,7,9 次谐波');subplot(2,2,4);plot(wt,y4,wt,y);title('1,3,5,19 次谐波');第2题参数如下：&g

45、t;> R=6e-4;C=17e-4;L=6e-9;模块参数设置：Integrator1 的 Initial condition:15kV在命令窗口为R,L,C赋值。仿真参数设置如下：Start time:0Stop time:100e-6Solver Type:Variable-step Solver:ode45Max step size:1e-7 Min step size:auto In itial step size:auto Relative tolera nce:1e-3 Absolute tolera nce:1e-6实验六 用matlab求解常微分方程1 .微分方程的概念

46、未知的函数以及它的某些阶的导数连同自变量都由一已知方程联系在一起的方程称为 微分方程。如果未知函数是一元函数，称为常微分方程。常微分方程的一般形式为F(t,y,y',y", ,y(n)0如果未知函数是多元函数，成为偏微分方程。联系一些未知函数的一组微分方程组称为微分方程组。微分方程中出现的未知函数的导数的最高阶解数称为微分方程的阶。若方程中未知函数及其各阶导数都是一次的，称为线性常微分方程，一般表示为y(n)ai(t)y(n1)an i(t)y' an (t)y b(t)若上式中的系数ai(t)，i1,2, ，n均与t无关，称之为常系数。2 常微分方程的解析解dy .

47、y 1有些微分方程可直接通过积分求解例如，一解常系数常微分方程dt可化为匹dtty 1,两边积分可得通解为 y ce 1.其中c为任意常数.有些常微分方程可用一些技巧,如分离变量法，积分因子法，常数变异法，降阶法等可化为可积分的方程而求得解析解.线性常微分方程的解满足叠加原理，从而他们的求解可归结为求一个特解和相应齐次微 分方程的通解.一阶变系数线性微分方程总可用这一思路求得显式解。高阶线性常系数微分 方程可用特征根法求得相应齐次微分方程的基本解，再用常数变异法求特解。一阶常微分方程与高阶微分方程可以互化，已给一个n阶方程(n)丄 /丄In(n 1) y f (t，y'，y"

48、，y )'(n 1)设y1y，* y， ，yn y，可将上式化为一阶方程组yr y2y2' yy'n 1ynyn' f(t，y1，y2，yn)反过来，在许多情况下，一阶微分方程组也可化为高阶方程。所以一阶微分方程组与高 阶常微分方程的理论与方法在许多方面是相通的，一阶常系数线性微分方程组也可用特征根法求解。3微分方程的数值解法除常系数线性微分方程可用特征根法求解，少数特殊方程可用初等积分法求解外，大部分微分方程无限世界，应用中主要依靠数值解法。考虑一阶常微分方程初值问题y'(t) f(t，y(t)，to t tfy(to) yo其中 y( yl, y2,

49、 ym )' , f ( f1, f2, fm)' , y0( yl0 , y20, ym0 )'-所谓数值解法，就是寻求y(t)在一系列离散节点totltntf上的近似值yk,k0,1,n称hktk 1tk为步长，通常取为常量h。最简单的数值解法是Euler法。Euler法的思路极其简单：在节点出用差商近似代替导数y'(tk)y(tki) y(tk)h这样导出计算公式(称为Euler格式)yk i yk hf (tk,yk), k 0,1,2,他能求解各种形式的微分方程。Euler法也称折线法。Euler方法只有一阶精度，改进方法有二阶Runge-Kutta法

50、、四阶Runge-Kutta法、五阶Runge-Kutta-Felhberg法和先行多步法等，这些方法可用于解高阶常微分方程(组)初值问题。边值问题采用不同方法，如差分法、有限元法等。数值算法的主要缺点是它缺乏物理理解。4 .解微分方程的 MATLAB令MATLAB主要用dsolve求符号解析解，ode45,ode23,ode15s求数值解。s=dsolve( 方程1', 方程2','初始条件1','初始条件2','自 变量)用字符串方程表示，自变量缺省值为t。导数用D表示，2阶导数用D2表示，以此类推。S返回解析解。在方程组情形，s为一

51、个符号结构。tout,yout=ode45( yprime ' ,t0,tf,y0)采用变步长四阶Runge-Kutta 法和五阶Runge-Kutta-Felhberg法求数值解，yprime是用以表示 f(t,y) 的M文件名，t0表示自变量的初始值，tf表示自变量的终值，y0表示初始向量值。输出向量tout表示节点(t 0,t 1,t n)T,输出矩阵yout表示数值解，每一列对应y的一个分量。若无输出参数，则自动作出图形。ode45是最常用的求解微分方程数值解的命令，对于刚性方程组不宜采用。ode23与ode45类似，只是精度低一些。ode12s用来求解刚性方程组，是用格式同o

52、de45。可以用helpdsolve, help ode45查阅有关这些命令的详细信息例1求下列微分方程的解析解(1)y'ay b(2)y''si n(2x)y,y(0)0,y'(0) 1(3)f'f g, g'g f,f'(0) 1, g'(0) 1方程(1)求解的MATLAB弋码为：>>clear;>>s=dsolve('Dy=a*y+b')结果为s =-b/a+exp(a*t)*C1方程(2)求解的MATLAB弋码为：>>clear;>>s=dsolve('D2y=si n(2*x)-y','y(0)=0','Dy(0)=1','x')>>simplify(s)%以最简形式显示s结果为s=(-1/6*cos(3*x)-1/2*cos(x)*si n(x)+(-1/2*si n(x)+1/6*si n(3*x)*cos(x)+5/3*si n(