中南大学系统仿真实验报告

1、本文格式为word版，下载可任意编辑中南大学系统仿真实验报告 ans = 试验一 matlab 中矩阵与多项式的基本运算 试验任务 1. 了解 matlab 命令窗口和程序文件的调用。 2 熟识如下 matlab 的基本运算： 矩阵的产生、数据的输入、相关元素的显示； 矩阵的加法、乘法、左除、右除； 特别矩阵：单位矩阵、 1 矩阵、 0 矩阵、对角阵、随机矩阵的产生和 运算； 多项式的运算：多项式求根、多项式之间的乘除。 基本命令训练 1、 eye(2) ans = 1 0 0 1 eye(4) ans = 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 2、 ones(2)

2、1 1 ans = 1 1 ones(4) ans = 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 ones(2,2) ans = 1 1 1 1 ones(2,3) ans = 1 1 1 1 1 1 ones(4,3) ans = 1 1 1 1 1 1 1 1 1 1 1 1 3、 zeros(2) 0 0 0 0 zeros(4) ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 zeros(2,2) ans = 0 0 0 0 zeros(2,3) ans = 0 0 0 0 0 0 zeros(3,2) ans = 0 0 0 0 00

3、 4、随机阵 rand(2,3) ans = 0.2785 0.9575 0.1576 0.5469 0.9649 0.9706 rand(2,3) ans = 0.9572 0.8003 0.4218 0.4854 0.1419 0.9157 5、 diag(5) ans = 5 diag(5,5) ans = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 5 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 diag(2,3) ans = 0 0 0 2 0 0 0 0 0 0 0 0 0 0 0 0 6、 (inv (a)为求 a 的逆矩阵)

4、b=5 3 1;2 3 8;1 1 1,inv(b) 5 3 1 2 3 8 1 1 1 ans = 0.6250 0.2500 -2.6250 -0.7500 -0.5000 4.7500 0.1250 0.2500 -1.1250 a=2 3;4 4,b=5 3;3 8,inv(a),inv(b);ab,a/b,inv(a)*b,b*inv(a) a = 2 3 4 4 b = 5 3 3 8 ans = -1.0000 0.7500 1.0000 -0.5000 ans = -2.7500 3.0000 3.5000 -1.0000 ans = 0.2258 0.2903 0.6452

5、0.2581 ans = -2.7500 3.0000 3.5000 -1.0000 ans = -2.0000 2.2500 5.0000 -1.7500 7、 p =1,-6,-72,-27, roots(p) p = 1 -6 -72 -27 ans = 12.1229 -5.7345 -0.3884 p=2,3,6,roots(p) p = 2 3 6 ans = -0.7500 + 1.5612i -0.7500 - 1.5612i 8、( a 为 n*n 的方阵) a=0 1 0;-4 4 0;-2 1 2,poly(a),b=sym(a),poly(b) a = 0 1 0 -4

6、 4 0 -2 1 2 ans = 1 -6 12 -8 b = 0, 1, 0 -4, 4, 0 -2, 1, 2 ans = x a 3-6*x a 2+12*x-8 9, 、( conv 是多项式相乘， deconv 是多项式相除) u=1 2 4 6 ,v=5 0 0 -6 7,conv(u,v) u = 1 2 4 6 v = 5 0 0 -6 7 ans = 5 10 20 24 -5 -10 -8 42 v=1 2 4 6 ,u=5 0 0 -6 7,deconv(u,v) v = 1 2 4 6 u = 5 0 0 -6 7 ans = 5 -10 10、（ 点乘是数组的运算，

7、没有点的乘是矩阵运算 ） a = 2 5;3 4, b =3 1;4 7,a.*b,a*b a = 2 5 3 4 b = 3 1 4 7 ans = 6 5 12 28 ans = 26 37 25 31 a = 2 3; b = 4 7; a.*b = 8 21; a*b %错误 a*b = 29; 11、（ who 可以看到你用过的一些变量， 来了） who your variables are: a b a ans b p u whos name size bytes whos 是把该变量及所存储的大小等信息都显示出 class attributes 2x2 32 double b 2

8、x2 32 double a 1x2 16 double ans 1x2 16 double b 1x2 16 double p 1x3 24 double u 1x5 40 double v 1x4 32 double 12、 a=2 5 3;6 5 4,disp(a),size(a),length(a) a = 2 5 3 6 5 4 2 5 3 6 5 4 ans = 2 3 ans = 3 试验二 matlab 绘图命令 试验任务 熟识 matlab 基本绘图命令，把握如下绘图方法： 1 坐标系的选择、图形的绘制； 2 图形注解（题目、标号、说明、分格线）的加入; 3 图形线型、符号、

9、颜色的选取 基本命令训练 1、t=0:pi/360:2*pi; x=cos(t)+ cos(t*4); y=si n( t)+ sin (t*4); xlabel(x 轴);ylabel(y 轴); plot(y,x),grid; 2、 t=0:0.1:100; x=3*t;y=4*t;z=si n( 2*t); plot3(x,y,z, g:) 15 i 0 5 0 05 1 1 5 2 3、x = linspace(-2*pi,2*pi,40); y=si n( x); stairs(x,y)4、 t=0:pi/360:2*pi; x=cos(t)+ cos(t*4) + sin (t*4

10、); y=si n( t)+ si n( t*4); plot(y,x, r:); xlabel(x 轴);ylabel(y 轴); 6、th=0:pi/20:2*pi; x=exp(j*th); plot(real(x),imag(x),r-.); grid; text(0,0,中心); 5、 th=0:pi/1000:2*pi; r=cos(2*th); polar(th,r); title( 四叶草图 ) 270 四叶草图 10 7、x=-2:0.01:2; 8、y=-2:0.01:2; 9、x,y = meshgrid(x,y); z = y.*exp(-x. a 2-y. a 2);

11、 c,h = con tour(x, yz); set(h,showtext,o n,textstep,get(h,levelstep)*2)_ 1os i ,5 2 n.s o o.s 8、x = 0:0.2:10; y = 2*x+3; subplot(411);plot(x,y); grid;title(y 的原函数); subplot(412) ;semilogy(x,y); grid;title(对 y 取对数); 丫的原画数 40 - 1 - 1 - 1 - 1 - 1 - 1 - 1 - l| i p | il | i| i 九 _ 1- _ _ i _ l _ : ： _ j

12、_ _ l _ u i| i |l i , il _ - j i i _ h- _ i i q 【 i 1 f i i i i ii i i it 10 1 1 2 3 4 5 6 r 6 9 10 , 对 y 取对数 对弋观对数 subplot(413) ;semilogx(x,y); 10 10 id 10 1lz 10 w 10 10 40 20 0 对好对数 2 把握循环、分支语句的编写，学会使用 look for 、 help 命令 grid;title(对 x 取对数); subplot(414) ;loglog(x,y);grid; title(对 xy 均取对数); 9、x =

13、 -3:0.3:3; bar(x,exp(-x.*x),g) 试验三 matlab 程序设计 试验任务 1 熟识 matlab 程序设计的方法和思路； 程序举例 1、 f=1,1; i=1; while f(i)+f(i+1)1000 f(i+2)=f(i)+f(i+1); i=i+1; end f,i f = columns 1 through 14 1 1 2 3 5 8 13 21 34 55 89 144 233 377 columns 15 through 16 610 987 i = 15 2、 m=3; n=4; for i=1:m for j=1:n a(i,j)=1/(i+j

14、-1); end end format rat 1 1/2 1/3 1/2 1/3 1/4 1/3 1/4 1/5 (分数格式形式。用有理数靠近显示数据 ) m=5; n=4; for i=1:m for j=1:n a(i,j)=1/(i+j-1); end end format rat a a = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/3 1/4 1/5 1/6 1/4 1/5 1/6 1/7 1/5 1/6 1/7 1/8 3、程序中没有 format rat 命令时，假如上次运行结果没有清除，输出的结果就是上次运行的 结果！但是运用 clear 命令清晰之前的

15、运行结果之后就会正常运行。 4、 x=input（请输入 x 的值:）; if x=10 y=cos(x+1)+sqrt(x*x+1); else y=x*sqrt(x+sqrt(x); 1/4 1/5 1/6 end y 请输入 x 的值 :2 y = 2391/647 x=input( 请输入 x 的值:); if x=10 y= fprintf ( 不在定义域内，请重新输入： );return else y=1/(x-10); end y 请输入 x 的值 :2 -1/8 5、 p=0 0 0 1 3 0 2 0 0 9; for i=1:length(p), if p(1)=0,p=p

16、(2:length(p); end ; end ; p p = columns 1 through 5 1 3 0 2 0 columns 6 through 7 0 9 p=0 0 0 1 3 0 2 0 0 9; p(p=0)=;p p = 1 3 2 9 6、 e2(500) ans = 1 1 2 3 5 8 13 21 34 55 89 144 233 377 lookfor ffibno e2 - ffibno 计算斐波那契亚数列的函数文件 help e2 ffibno 计算斐波那契亚数列的函数文件 n 可取任意自然数 程序如下 （用法： lookfor 关键词 在全部 m 文件中

17、找 关键词 ，比如： lookfor max （即查找关键词 max） 其实就和我们平常用 ctrl+f 来查找关键词是一样的 而 help 是显示 matlab 内置的关心信息 用法： help 命令，比如 help inv ，作用就是调用 inv 这个命令的关心） 程序设计题 用一个 matlab 语言编写一个程序：输入一个自然数，推断它是否是素数， 假如是，输出 it is one prime ，假如不是，输出 it is n ot o ne prime. 。要求通 过调用子函数实现。 最好能具有如下功能： 设计较好的人机对话界面， 程序中 含有提示 性的输入输出语句。 能实现循环操作，

18、 由操输入相关命令来掌握 是否连续进行素数的推断。 假如操盼望停止这种推断， 则可以退出程序。 假如所输入的自然数是一个合数， 除了给出其不是素数的结论外， 还应给出至少 一种其因数分解形式。例：输入 6 ，由于 6 不是素数。则程序中除了有 it is not one prime 的结论外，还应有： 6=2*3 的说明。 function sushu while 1 x=input（ 请输入一个自然数 ）; if x2 disp（ 既不是质数又不是合数 ）; else if isprime(x)=1 disp（ 这是一个素数 ）; else disp( 这是一个合数，可以因式分解为： ) f

19、or n=2:sqrt(x) if rem(x,n)=0 num3=x; num1=n; num2=x/n; disp(num2str(num3),=,num2str(num1),x,num2str(num2) end end end end y=input( 是否连续推断？连续请按 1 ,按任意键退出 ： ) ； if y=1 break end end 试验四 matlab 的符号计算与 simulink 的使用 试验任务 1 .把握 matlab 符号计算的特点和常用基本命令; 2 .把握 simulink 的使用。 程序举例 1 . 求矩阵对应的行列式和特征根 a=sym(a11 a1

20、2;a21 a22); da=det(a) ea=eig(a) da = a11*a22-a12*a21 ea = 1/2*a11+1/2*a22+1/2*(a11 a 2-2*a11*a22+a22 a 2+4*a12*a21) a (1 /2) 1/2*a11+1/2*a22-1/2*(a11 a 2-2*a11*a22+a22 a 2+4*a12*a21) a (1 /2) a=sym(2 3;1 5); da=det(a) ea=eig(a) da = 7 ea = 7/2+1/2*21a(1/2) 7/2-1/2*21a(1/2) 2. 求方程的解(包括精确解和肯定精度的解) r1=

21、solve(xa2+x-1) rv=vpa(r1) rv4=vpa(r1,4) rv20=vpa(r1,20) r1 = 1/2*5 八( 1/2)-1/2 -1/2*5 八( 1/2)-1/2 rv = .6894848260 -1.689484826 rv4 = .6180 -1.618 rv20 = .689484820 -1.68948482 3 a=sym( a );b=sym( b );c=sym( w=10;x=5;y=-8;z=11; a=a,b;c,d b=w,x;y,z det(a) det(b) a = a, b c );d=sym( d ); %定义 4 个符号变量 %

22、定义 4 个数值变量 %建立符号矩阵 a %建立数值矩阵 b %计算符号矩阵 a 的行列式 %计算数值矩阵 b 的行列式 c, d x 10 -8 11 ans = a*d-b*c ans = 150 4. syms x y;s=(-7*x a 2-8*y a 2)*(-x a 2+3*y a 2); expand(s) % 对 s 绽开 collect(s,x) % 对 s 按变量 x 合并同类项 （ 无同类项 ） factor(ans) % 对 ans 分解因式 ans = 7*xa4-13*xa2*ya2-24*ya4 ans = 7*xa4-13*xa2*ya2-24*ya4 ans

23、= (8*ya2+7*xa2)*(xa2-3*ya2) 5. 对方程 ax=b 求解 a=34,8,4;3,34,3;3,6,8; b=4;6;2; x=linsolve(a,b) % 调用 linsolve 函数求解 ab % 用另一种方法求解 0.0675 0.1614 diff(f) %未指定求导变量和阶数，按缺省规章处理 0.1037 ans = 0.0675 0.1614 0.1037 6 对方程组求解 a11*x1+a12*x2+a13*x3=b1 a21*x1+a22*x2+a23*x3=b2 a31*x1+a32*x2+a33*x3=b3 syms a11 a12 a13 a2

24、1 a22 a23 a31 a32 a33 b1 b2 b3; a=a11,a12,a13;a21,a22,a23;a31,a32,a33; b=b1;b2;b3; xx=ab %用左除运算求解 (x=linsolve(a,b) %调用 linsolve 函数求的解 ) xx = (a12*a23*b3-a12*b2*a33+a13*a32*b2-a13*a22*b3+b1*a22*a33-b1*a32*a23)/( a11*a22*a33-a11*a32*a23-a12*a21*a33+a32*a21*a13-a22*a31*a13+a31*a12*a 23) -(a11*a23*b3-a1

25、1*b2*a33-a21*a13*b3-a23*a31*b1+b2*a31*a13+a21*b1*a33)/ (a11*a22*a33-a11*a32*a23-a12*a21*a33+a32*a21*a13-a22*a31*a13+a31*a12* a23) (a32*a21*b1-a11*a32*b2+a11*a22*b3-a22*a31*b1-a12*a21*b3+a31*a12*b2)/( a11*a22*a33-a11*a32*a23-a12*a21*a33+a32*a21*a13-a22*a31*a13+a31*a12*a 23) 7 syms a b t x y z; f=sqrt

26、(1+exp(x); %求 f 对 x 的二阶导数 %求 f 对 x 的三阶导数 %按参数方程求导公式求 y 对 x 的导数 ans = 1/2/(1+exp(x) a (1/2)*exp(x) ans = -2*si n(x)-x*cos(x) ans = -3*cos(x)+x*si n(x) ans = -b*cos(t)/a/si n(t) 三、 simulink 的使用 f=x*cos(x); diff(f,x,2) diff(f,x,3) f1=a*cos(t);f2=b*si n(t); diff(f2)/diff(f1) diff(f) %未指定求导变量和阶数，按缺省规章处理

27、仿真图: 波形图: 其中: g i ( r 0.14 试验五 matlab 在掌握系统分析中的应用 试验任务 1. 把握 matlab 在掌握系统时间响应分析中的应用; 2. 把握 matlab 在系统根轨迹分析中的应用； 3. 把握 matlab 掌握系统频率分析中的应用； 4. 把握 matlab 在掌握系统稳定性分析中的应用 基本命令 1. step 2. impulse 3. in itial 4. isim 5. rlocfi nd 6. bode 7. margin 8. nyquist 9. nichols 10. cloop 程序举例 1. 求下面系统的单位阶跃响应 1.5 g

28、(s) t num=4 ; den=1 , 1 , 4; step( num , den) y , x , t=step( num , den); tp=spli ne(y , t , max(y) % 计算峰值时间 max(y) % 计算峰值 tp = step response 0.5 time (sec) p m 1.6062 ans = 1.4441 0.18 0.16 2. 求如下系统的单位阶跃响应 x i 0 1 x i 0 u x 2 6 5 x 2 1 a=0,1;-6,-5;b=0;1;c=1,0;d=0; y,x=step(a,b,c,d); plot(y) 1,0 x i

29、 x 2 3. 求下面系统的单位脉冲响应: g(s) 儿 num=4 ; den=1 , 1 ,4; impulse( nu m,de n) time (sec) response to initial conditions 4. 已知二阶系统的状态方程为: e u p m a 0.25 0.2 0.15 0.1 0.05 0 -0.05 l i 1 - - r - 1 厂、 / impulse response c=1 , 0 ； d=0; xo=1 ,0 ； subplot(1 , 2,1); in itial(a , b , c ,d,x0) subplot(1 , 2,2); impu

30、lse(a , b , c , d) 5 :系统传递函数为: g(s) 占 输入正弦信号时，观看输出 信号的相位差。 num=1 ; den=1 ,1; t=0 : 0.01 : 10 ; u=s in( 2*t) ; hold on plot(t,u, r) lsim( nu m,de n,u,t) real axis 6. 有一二阶系统，求出周期为 4 秒的方波的输出响应 2s 2 5s 1 2 s 2s 3 num=2 5 1; den=1 2 3; t=(0:.1:10); period=4; u=(rem(t,period)=period./2); % 看 rem 函数功能 lsim

31、( nu m,de n,u,t); 7. 已知开环系统传递函数，绘制系统的根轨迹，并分析其稳定性 g(s) k(s 2) (s 2 4s 3) 2 num=1 2; den 仁 1 4 3; den=conv (de n1,de n1); figure(1) -6 rlocus(num,den) k,p= rlocfi nd(nu m,de n) l l l 匚一 l c - - - r t r r -10 -8 -6 -4 -2 4 6 root locus 0 2 2 2 0 0 2 2 - - cpxa vyanma n -8 2.5 g(s) lin ear simulati on r

32、esults1.5 e 0.5 -0.5 -1 -1.5 -2 time (sec) p m impulse response (k=55) 5 5 figure(2) k=55; num 仁 k*1 2; den=1 4 3; den 1=c onv (de n,den); nu m,de n=cloop( nu m1,de n1,-1); impulse( nu m,de n) title(impulse resp onse (k=55) -1.5 -1 0 200 400 600 time (sec) 800 1000 1200 5 5 0 0 5 5 0 0 o - - 6 x 10

33、8 ! - impulse resp on se(k=56) figure(3) k=56; num 仁 k*1 2; den=1 4 3; den 1=c onv (de n,den); nu m,de n=cloop( nu m1,de n1,-1); impulse( nu m,de n) -2 -4 -6 -8 l 0 500 1000 time (sec) title(impulse resp on se(k=56) 1500 2021 2500 select a point in the graphics win dow selected_po int = -2.5924 - 0.

34、0248i 0.7133 -3.4160 -2.5918 -0.9961 + 0.4306i -0.9961 - 0.4306i bode diagram 8. 作如下系统的 bode 图 g(s) n=1 , 1 ; d=1 , 4,11 , 7; bode( n , d),grid on frequency (rad/sec) s 1 s 3 4s 2 11s 7 9. 系统传函如下 g(s) s 1 0.5s e (s 2) 3 求有理传函的频率响应，然后在同一张图上绘出以四阶伯德近似表示的系统频 率响应 num=1;de n=co nv(1 2,co nv(1 2,1 2); -20

35、w=logspace(-1,2); t=0.5; ylabel(gai n); subplot(2,1,2); semilogx(w,p1,w,p2,g-);grid on; xlabel(freque ncy); ylabel(phase); 10. 已知系统模型为 求它的幅值裕度和相角裕度 n=3.5; d=1 2 3 2; gm,pm,wcg,wcp=margi n(n,d) g(s) 3.5 s 3 2s 2 3s 2 bode plot m1,p1=bode( nu m,de n, 2); p1=p1-t*w*180/pi; n 2,d2=pade(t,4); nu mt=c onv

36、(n2,nu m); den t=(c onv (de n, d2); 2 freque xlabel(freque gm = 1.1433 pm = 7.1688 wcg = 1.7323 wcp = 1.6541 nyq uist (n, d1); hold on nyq uist (n, d2) ; nyq uist (n, d3) ; nyq uist (n, d4); nyquist diagram 80 | - - - - 11. 二阶系统为: g(s) n s 2 2 n s s 令 wn=1, 分别作出 e =2,1 , 0.707 , c a 5 y a 时的 nyquist

37、 曲线。 m n=1; d 仁 1 , 4,1 ; d2=1 , 2 , 1 ； d3=1 , 1.414,1; d4=1,1,1; nyquist diagram real axis 12. 已知系统的开环传递函数为 -20 _一 - 一 一 -一_一_ - 60 00 5 5 s s 2) s s 3 3 2 2 s s v vr ra an nk ky ya a 14. 一多环系统，其结构图如下，使用 nyquist 频率曲线推断系统的稳定性。 16.7s (0.85s 1)(0.25s 1)(0.0625s 1)绘制系统的 nyqusit 图, 并争论系统的稳定性 . g=tf(100

38、0,co nv(1,3,2,1,5); nyquist(g);axis(square) 13. 分别由 w 的自动变量和人工变量作下列系统的 nyquisf 曲线： m g(s) 1 s(s 1) n=1 ; d=1 , 1 ,0; nyquist(n ,d) ; % 自动变量 n=1 ; d=1 , 1 ,0; w=0.5 : 0.1 : 3; nyq uist (n , d , w) ; % 人工变量 -2 -1 -1.5 -0.9 -0.8 -0.7 -0.6 -0.5 real axis -0.4 -0.3 -0.2 -0.1 5 5 - -0 0 s si ix xa a u u a

39、 an nt t9 9a a 卩 g(s) 40 30 20 10 0 figure(2) nu m2,de n2=cloop( nu m,de n); impulse( nu m2,de n2); k1=16.70.0125;z1=0; p1=-1.25 -4 -16; nu m1,de n1=zp2tf(z1,p1,k1); nu m,de n=cloop( nu m1,de n1); z,p,k=tf2zp( nu m,de n);p figure(1) nyq uist (nu m,de n) -2 -1 nyquist diagram -0.5 0 0.5 real axis 1 1

40、.5 1 0.5 s a 0 n g m -0.5 -1 -1.5 1.5 20 impulse resp onse -10.5969 +36.2148i -10.5969 -36.2148i -0.0562 15. 已知系统为: eanmlpm -5 -10 -15 0 0.1 0.2 0.3 time (sec) nichols chart 0.4 0.5 0.6 10 -2 10 0 10 2 open-loop phase (deg) frequency (rad/sec) g(s) s(s 1) 作该系统的 nichols 曲线。 n=1 ; d=1 , 1 , 0; ni chol

41、s( n , d); 16. 已知系统的开环传递函数为: g(s) k s(s 1)(s 2) 当 k=2 时，分别作 nichols 曲线和波特图 num=1; den=conv (co nv(1 0,1 1),0.5 1); subplot(1,2,1); ni chols( nu m,de n);grid; b % n ichols 曲线 g subplot(1,2,2); 2 g=tf(nu m,de n); bode(feedback(g,1,-1);grid; nichols chart koflr ea mhnaam msuvesab 50 50 - - 00 80 -270 b

42、ode diagram 90 - - % 波特图 17. 系统的开环传递函数为: 分别确定 k=2 和 k=10 时闭环系统的稳定性 d 仁 1 , 3,2,0 ; n 仁 2; nc1 , dc1=cloop( n1 , d1 ,-1); roots(dc1) d2=d1 ; n2=10; nc2 , dc2=cloop(n2 , d2,-1) ; roots(dc2) ans = -2.5214 -0.2393 + 0.8579i -0.2393 - 0.8579i ans = -3.3089 0.1545 + 1.7316i 0.1545 - 1.7316i 18. 系统的状态方程为：

43、x 1 4 3 0 x 1 1 x 2 1 0 0 x 2 0 u x 3 0 1 0 x 3 0 x 1 y 0 1 2 x 2 x 3 g(s) k s(s 1)(s 2) 试确定系统的稳定性。 a=-4,-3,0 ; 1,0,0 ; 0,1,0 ; b=1;0;0 ; c=0,1,2 ; d=0 ; eig(a) % 求特征根 ran k(ctrb(a,b) ans = 0 -1 -3 ans = 3 试验六连续系统数字仿真的基本算法 试验任务 1 .理解欧拉法和龙格 - 库塔法的基本思想； 2 理解数值积分算法的计算精度、速度、稳定性与步长的关系; 程序举例 1. 取 h=0.2 ,试

44、分别用欧拉法、 rk2 法和 rk4 法求解微分方程的数值解，并 比较计算精度。 注：解析解：y(t) .1 2t clear t(1)=0 ； y(1)=1; y_euler(1)=1; y_rk2(1)=1; y_rk4(1)=1; h=0.001; % 步长修改为 0.001 for k=1:5 y(t) y(t) y(0) 2t y(t) 1 t(k+1)=t(k)+h; y(k+1)=sqrt(1+2*t(k+1); end for k=1:5 y_euler(k+1)=y_euler(k)+h*(y_euler(k)-2*t(k)/y_euler(k); end for k=1:5

45、 k1= y_rk2(k)-2*t(k)/y_rk2(k); k2=(y_rk2(k)+h*k1)-2*(t(k)+h)/(y_rk2(k)+h*k1); y_rk2(k+1)=y_rk2(k)+h*(k1+k2)/2; end for k=1:5 k1= y_rk4(k)-2*t(k)/y_rk4(k); k2=(y_rk4(k)+h*k1/2)-2*(t(k)+h/2)/(y_rk4(k)+h*k1/2); k3=(y_rk4(k)+h*k22)-2*(t(k)+h/2)/(y_rk4(k)+h*k2/2); k4=(y_rk4(k)+h*k3)-2*(t(k)+h)/(y_rk4(k)+

46、h*k3); y_rk4(k+1)=y_rk4(k)+h*(k1+2*k2+2*k3+k4)/6; end disp( 时间 解析解 yt=t, y, y_euler, y_rk2, y_rk4; disp(yt) 欧拉法 rk2 法 rk4 法 ) 时间 解析解 欧拉法 rk2 法 rk4法0 1.0000 1.0000 1.0000 1.0000 0.0010 1.0010 1.0010 1.0010 1.0010 0.0020 1.0020 1.0020 1.0020 1.0020 0.0030 1.0030 1.0030 1.0030 1.0030 0.0040 1.0040 1.00

47、40 1.0040 1.0040 0.0050 1.0050 1.0050 1.0050 1.0050 y(t) 2ry(t) y(t) 0 在 0 t 10 上的数字仿真解(已知： y(0) y(0) 0)， 并将不同步长下的仿真结果与解析解进行精度比较。 说明： 已知该微分方程的解析解分别为： 100， yt y t 100cost ( 当 r 0) 100e 2t cs 仝 t 10 e*s in 仝 t 2 3 2 ( 当 r 0.5) 采纳 rk4 法进行计算，选择状态变量: 2. 考虑如下二阶系统： x 1 y x 2 y 则有如下状态空间模型及初值条件 x 1 x 2 x 1 (0) 100 x 2 x 1 2rx 2 x 2 (0) 0 y x 1 采纳 rk4 法进行计算。 clear h=input( 请输入步长 h=); m=round(10/h); t(1)=0; y_0(1)=100; y_05(1)=100; 和 y_05 分别对应于为 r=0 和 r=0.5) x1(1)=100; x2(1)=0; y_rk4_0(1)=x