MATLAB刘卫国版本

海 南 大 学

MATLAB 实验 (刘卫国版本)

题 目： MATLAB 实验 学 号： 2009050******

姓 名： **** 年 级： 09电气 学 院： 机电工程学院 系 别： 电气工程系 专 业： 电气工程及其自动化

1

实验一 MATLAB运算基础

二（1）二、实验内容

1.先求下列表达式的值，然后显示MATLAB工作空间的使用情况并保存全部变量。 ⑴z2sin85?1?1?e2；

程序为>> z1=2*sin(85*pi/180)/(1+exp(2)) z1 =

0.2375 ⑵z12ln(x?1?x2)，其中x???21?2i?2???0.455??； >> x=[2 1+2i;-0.45 5];

>> z2=1/2*log(x+sqrt(1+x^2)) z2 =

0.7114 - 0.0253i 0.8968 + 0.3658i 0.2139 + 0.9343i 1.1541 - 0.0044i

⑶ze0.3a?e?0.3a3?2sin(a?0.3)?ln0.3?a2,a??3.0,?2.9,?2.8,?,2.8,2.9,3.0程序为>> a=(-3.0:0.1:3.0);

>> z3=(exp(0.3.*a)-exp(-0.3.*a))./2.*sin(a+0.3)+log((0.3+a)./2) z3 =

Columns 1 through 3

0.7388 + 3.1416i 0.7696 + 3.1416i 0.7871 + 3.1416i

Columns 4 through 6

0.7913 + 3.1416i 0.7822 + 3.1416i 0.7602 + 3.1416i

2

Columns 7 through 9

0.7254 + 3.1416i 0.6784 + 3.1416i 0.6196 + 3.1416i

Columns 10 through 12

0.5496 + 3.1416i 0.4688 + 3.1416i 0.3780 + 3.1416i

Columns 13 through 15

0.2775 + 3.1416i 0.1680 + 3.1416i 0.0497 + 3.1416i

Columns 16 through 18

-0.0771 + 3.1416i -0.2124 + 3.1416i -0.3566 + 3.1416i

Columns 19 through 21

-0.5104 + 3.1416i -0.6752 + 3.1416i -0.8536 + 3.1416i

Columns 22 through 24

-1.0497 + 3.1416i -1.2701 + 3.1416i -1.5271 + 3.1416i

Columns 25 through 27

-1.8436 + 3.1416i -2.2727 + 3.1416i -2.9837 + 3.1416i

Columns 28 through 30

-37.0245 -3.0017 -2.3085

Columns 31 through 33

-1.8971 -1.5978 -1.3575

Columns 34 through 36

-1.1531 -0.9723 -0.8083

Columns 37 through 39

-0.6567 -0.5151 -0.3819

3

Columns 40 through 42

-0.2561 -0.1374 -0.0255

Columns 43 through 45

0.0792 0.1766 0.2663

Columns 46 through 48

0.3478 0.4206 0.4841

Columns 49 through 51

0.5379 0.5815 0.6145

Columns 52 through 54

0.6366 0.6474 0.6470

Columns 55 through 57

0.6351 0.6119 0.5777

Columns 58 through 60

0.5327 0.4774 0.4126

Column 61

0.3388

?t2 ⑷z?0?t?1?t24??1 1?t?2，其中t?0:0.5:2.5

??t2?2t?1,2?t?3程序为>> t=(0:0.5:2.5);

>> z4=(t>=0&t<1).*(t.^2)+(t>=1&t<2).*(t.^2-1)+(2<=t&t<3).*(t.^2-2.*t+1) z4 =

0 0.2500 0 1.2500 1.0000 2.2500

4

?1234?4??13?1?????，求下列表达式的值：

32.已知：A?34787,B?20???????3657???3?27?? ⑴A?6*B和A?B?I(其中I为单位矩阵)；⑵A*B和A.*B; ⑶A^3和A.^3;⑷A/B和A\\B;⑸[A,B]和[A([1,3]),:);B^2]. 程序为>> A=[12 34 -4;34 7 87;3 65 7];

>> B=[1 3 -1;2 0 3;3 -2 7]; >> I=[1 0 0;0 1 0;0 0 1]; >> A+6*B

ans =

18 52 -10 46 7 105 21 53 49

>> A-B+I

ans =

12 31 -3 32 8 84 0 67 1

>> A^3

ans =

37226 233824 48604 247370 149188 600766 78688 454142 118820

>> A.^3

ans =

1728 39304 -64 39304 343 658503 27 274625 343

>> A/B

5

ans =

16.4000 -13.6000 7.6000 35.8000 -76.2000 50.2000 67.0000 -134.0000 68.0000

>> B\\A

ans =

109.4000 -131.2000 322.8000 -53.0000 85.0000 -171.0000 -61.6000 89.8000 -186.2000

>> [A,B]

ans =

12 34 -4 1 3 -1 34 7 87 2 0 3 3 65 7 3 -2 7

>> [A([1,3],:);B^2]

ans =

12 34 -4 3 65 7 4 5 1 11 0 19 20 -5 40

?12345??30?678910??17?6???3.设有矩阵A和B：A??1112131415?,B??023???1617181920???97???2122232425???413⑴求它们的乘积C；⑵将矩阵C的右下角3?2子矩阵赋给D；

16?9???4? ?0?11??⑶查看MATLAB工作空间的使用情况。

程序为>> 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 25]; >> B=[3 0 16;17 -6 9 ;0 23 -4;9 7 0;4 13 11]; >> C=A*B

C = 93 150 77

6

258 335 237 423 520 397 588 705 557 753 890 717 >> D=C([3,4,5],[2,3])

D = 520 397 705 557 890 717 >> whos

Name Size Bytes Class Attributes A 5x5 200 double B 5x3 120 double C 5x3 120 double D 3x2 48 double 4.完成下列操作：

⑴求[100,999]之间能被21整除的数的个数； ⑵建立一个字符串向量，删除其中的大写字母。 程序为>> A=100:1:999; K=find(rem(A,21)==0); length(K) ans = 43

>> ch='aHDJLK143663CFHI'; >> H=find(ch>='A'&ch<'Z'); >> ch(H)=[ ] ch = a143663

实验二 MATLAB矩阵分析与处理二（1）E=eye(3);

R=rand(3,2); O=zeros(2,3); S=diag([2,3]); A=[E,R;O,S]; A^2

B=[E,(R+R*S);O,S^2] ans =

1.0000 0 0 2.8504 1.9439

7

0 1.0000 0 0.6934 3.5652 0 0 1.0000 1.8205 3.0484 0 0 0 4.0000 0 0 0 0 0 9.0000 B =

1.0000 0 0 2.8504 1.9439 0 1.0000 0 0.6934 3.5652 0 0 1.0000 1.8205 3.0484 0 0 0 4.0000 0 0 0 0 0 9.0000

（2）>> H=hilb(5) P=pascal(5) Hh=det(H) Hp=det(P) Th=cond(H) Tp=cond(P) H =

1.0000 0.5000 0.3333 0.2500 0.2000 0.5000 0.3333 0.2500 0.2000 0.1667 0.3333 0.2500 0.2000 0.1667 0.1429 0.2500 0.2000 0.1667 0.1429 0.1250 0.2000 0.1667 0.1429 0.1250 0.1111 P =

1 1 1 1 1 1 2 3 4 5 1 3 6 10 15 1 4 10 20 35 1 5 15 35 70 Hh =

3.7493e-012

8

Hp = 1 Th =

4.7661e+005 Tp =

8.5175e+003

（3）>> A=fix(10*rand(5)) H=det(A) Trace=trace(A) Rank=rank(A) Norm=norm(A) A =

7 8 0 7 6 4 5 6 4 3 9 2 3 3 5 4 6 8 1 1 4 8 5 1 6 H =

-6498

Trace = 22 Rank = 5

9

Norm =

23.8478

（4）>> A=[-29,6,18;20,5,12;-8,8,5] [V,D]=eig(A) %数学意义略 A =

-29 6 18 20 5 12 -8 8 5 V =

0.7130 0.2803 0.2733 -0.6084 -0.7867 0.8725 0.3487 0.5501 0.4050 D =

-25.3169 0 0 0 -10.5182 0 0 0 16.8351

（5）A=hilb(4) A(:,1)=[] A(4,:)=[]

B=[0.95,0.67,0.52]'; X=inv(A)*B

B1=[0.95,0.67,0.53]'; X1=inv(A)*B1 N=cond(B) N1=cond(B1)

Na=cond(A) %矩阵A为病态矩阵 A =

1.0000 0.5000 0.3333 0.2500

10

0.5000 0.3333 0.2500 0.2000 0.3333 0.2500 0.2000 0.1667 0.2500 0.2000 0.1667 0.1429 A =

0.5000 0.3333 0.2500 0.3333 0.2500 0.2000 0.2500 0.2000 0.1667 0.2000 0.1667 0.1429 A =

0.5000 0.3333 0.2500 0.3333 0.2500 0.2000 0.2500 0.2000 0.1667 X =

1.2000 0.6000 0.6000 X1 =

3.0000 -6.6000 6.6000 N = 1 N1 = 1

11

Na =

1.3533e+003

（6）A=[1,4,9;16,25,36;49,64,81] B=sqrtm(A)

C=sqrt(A) %sqrtm函数是以矩阵为单位进行计算，sqrt函数是以矩阵中的元素进行计算 A =

1 4 9 16 25 36 49 64 81 B =

0.6344 + 1.3620i 0.3688 + 0.7235i 0.7983 - 0.4388i 1.4489 + 1.1717i 2.7697 + 0.6224i 3.2141 - 0.3775i 4.3578 - 1.6237i 5.7110 - 0.8625i 7.7767 + 0.5231i C =

1 2 3 4 5 6 7 8 9

实验三 选择程序结构设计

二、（1） x=[-5.0,-3.0,1.0,2.0,2.5,3.0,5.0]; y=[]; %建立存放所有y值的矩阵 for x0=x

if x0<0&x0~=-3 y=[y,x0*x0+x0-6];

elseif x0>=0&x0<5&x0~=2&x0~=3 y=[y,x0*x0-5*x0+6]; else

y=[y,x0*x0-x0-1]; end end

x %输出所有x y

12

x =

-5.0000 -3.0000 1.0000 2.0000 2.5000 3.0000 5.0000 y =

14.0000 11.0000 2.0000 1.0000 -0.2500 5.0000 19.0000

（2）x=input('请输入一个百分制成绩：'); if x>100|x<0

disp('您输入的成绩不是百分制成绩，请重新输入。'); else

if x<=100&x>=90 disp('A');

elseif x<=89&x>=80 disp('B');

elseif x<=79&x>=70 disp('C');

elseif x<=69&x>60 disp('D'); else disp('E'); end end

请输入一个百分制成绩：89 B

（3）n=input('请输入员工工号：'); h=input('该员工工作时数是：'); if h>120

x=(h-120)*84*(1+0.15)+120*84; elseif h<60 x=h*84-700; else x=h*84; end

disp([num2str(n),'号员工','的应发工资为',num2str(x)]); 请输入员工工号：59 该员工工作时数是：45 59号员工的应发工资为3080

13

（4）>> a=fix(10+(99-10)*rand(1,2)) %产生两个随机整数 x=a(1); y=a(2);

t=input('请输入运算符号：','s'); if t=='+' z=x+y; elseif t=='-' z=x-y; elseif t=='*' z=x*y; elseif t=='/' z=x/y; end

disp([num2str(x),t,num2str(y),'=',num2str(z)]) a =

43 86

请输入运算符号：+ 43+86=129

（5）>> a=rand(5,6) %产生5x6的随机矩阵 n=input('请输入您要输出矩阵的第几行：'); if n>5

disp('超出了矩阵的行数，矩阵的最后一行为：') a(5,:) else

disp(['矩阵的第',num2str(n),'行为：']) a(n,:) end a =

0.8537 0.6449 0.3412 0.5681 0.6946 0.8801 0.5936 0.8180 0.5341 0.3704 0.6213 0.1730 0.4966 0.6602 0.7271 0.7027 0.7948 0.9797 0.8998 0.3420 0.3093 0.5466 0.9568 0.2714 0.8216 0.2897 0.8385 0.4449 0.5226 0.2523

14

请输入您要输出矩阵的第几行：85 超出了矩阵的行数，矩阵的最后一行为： ans =

0.8216 0.2897 0.8385 0.4449 0.5226 0.2523

实验四

二（1）>> s=0; n=input('n=?'); for i=1:n s=s+1/i/i; end

PI=sqrt(6*s) pi n=? PI = 0 ans =

3.1416

（2）>> y=0; n=1; while(y<3) y=y+1/(2*n-1); n=n+1; end

y=y-1/(2*(n-1)-1) n=n-2 y =

2.9944

循环结构程序设计

15

n =

56

（3）>> a=input('a=?'); b=input('b=?'); Xn=1;

Xn1=a/(b+Xn); n=0;

while abs(Xn1-Xn)>1e-5 Xn=Xn1; Xn1=a/(b+Xn); n=n+1; if n==500 break; end end n Xn1

r1=(-b+sqrt(b*b+4*a))/2 r2=(-b-sqrt(b*b+4*a))/2 a=?85 b=?56 n = 3 Xn1 =

1.4788 r1 =

1.4788 r2 =

16

-57.4788

（4）for i=1:100 if i==1 f(i)=1; elseif i==2 f(i)=0; elseif i==3 f(i)=1; else

f(i)=f(i-1)-2*f(i-2)+f(i-3); end end max(f) min(f) sum(f)

length(find(f>0)) length(find(f==0)) length(find(f<0)) ans =

4.3776e+011 ans =

-8.9941e+011 ans =

-7.4275e+011 ans = 49 ans = 2

17

ans = 49

（5）>> s=0;n=0; for i=2:49 b=i*(i+1)-1; m=fix(sqrt(b)); for j=2:m if rem(b,j)==0 break end end if j==m n=n+1; s=s+b; end end n s n = 28 s =

21066

实验五 函数文件

二（1）函数fushu.M文件：

function [e,l,s,c] = fushu(z)

%fushu 复数的指数，对数，正弦，余弦的计算 %e 复数的指数函数值 %l 复数的对数函数值 %s 复数的正弦函数值

18

%c 复数的余弦函数值 e=exp(z); l=log(z); s=sin(z); c=cos(z);

z=input('请输入一个复数z='); [a,b,c,d]=fushu(z) 请输入一个复数z=1+i 命令文件M：

z=input('请输入一个复数z='); [a,b,c,d]=fushu(z) a =

1.4687 + 2.2874i b =

0.3466 + 0.7854i c =

1.2985 + 0.6350i d =

0.8337 - 0.9889i

19

（2）m1=input('输入m1='); m2=input('输入m2='); theta=input('输入theta='); x=theta*pi/180; g=9.8;

A=[m1*cos(x) -m1 -sin(x) 0 m1*sin(x) 0 cos(x) 0 0 m2 -sin(x) 0 0 0 -cos(x) 1]; B=[0;m1*g;0;m2*g]; X= function (A,B) 输入m1=1 输入m2=1 输入theta=30 X =

7.8400 3.3948 6.7896 15.6800

（3）function flag=mat3(x) flag=1;

for i=2:sqrt(x) if rem(x,i)==0 flag=0; break; end end

%在命令窗口调用该函数文件： for i=10:99

j=10*rem(i,10)+fix(i/10); if mat3(i)&mat3(j)

20

disp(i) end end j =

11 31 71 13 73 17 37 97 79

（4）函数fx.m文件： function f= fx(x)

%fx fx求算x矩阵下的f(x)的函数值 A=0.1+(x-2).^2; B=0.01+(x-3).^4; f=1./A+1./B;

命令文件： clc;

x=input('输入矩阵x='); f=fx(x)

>> x=input('输入矩阵x='); f=fx(x) 运算结果

输入矩阵x=[7 2;12 5] f =

0.0437 10.9901 0.0101 0.1724

（5）（1）函数f.m文件: function f=f(x) f=x+10*log(x^2+5);

命令文件： clc;

n1=input('n1='); n2=input('n2='); n3=input('n3='); y1=f(n1); y2=f(n2); y3=f(n3);

21

y=y1/(y2+y3) 运算结果如下 n1=40 n2=30 n3=20 y =

0.6390

（5）（2）函数g.m文件 function s= g(n) for i=1:n g(i)=i*(i+1); end s=sum(g);

命令文件： clc;

n1=input('n1='); n2=input('n2='); n3=input('n3='); y1=g(n1); y2=g(n2); y3=g(n3); y=y1/(y2+y3) 运算结果如下 n1=40 n2=30 n3=20 y =

1.7662

实验六 高层绘图操作（1）>> x=linspace(0,2*pi,101);

y=(0.5+3*sin(x)./(1+x.^2)).*cos(x); plot(x,y)

22

（2）a >> x=linspace(-2*pi,2*pi,100); y1=x.^2; y2=cos(2*x); y3=y1.*y2;

plot(x,y1,'b-',x,y2,'r:',x,y3,'y--'); text(4,16,'\\leftarrow y1=x^2');

text(6*pi/4,-1,'\\downarrow y2=cos(2*x)'); text(-1.5*pi,-2.25*pi*pi,'\%uparrow y3=y1*y2');

23

B x=linspace(-2*pi,2*pi,100); y1=x.^2; y2=cos(2*x); y3=y1.*y2;

subplot(1,3,1);%分区 plot(x,y1);

title('y1=x^2');%设置标题 subplot(1,3,2); plot(x,y2);

title('y2=cos(2*x)'); subplot(1,3,3); plot(x,y3);

title('y3=x^2*cos(2*x)');

24

C >> x=linspace(-2*pi,2*pi,20); y1=x.^2;

subplot(2,2,1);%分区 bar(x,y1);

title('y1=x^2的条形图');%设置标题 subplot(2,2,2); stairs(x,y1);

title('y1=x^2的阶梯图'); subplot(2,2,3); stem(x,y1);

title('y1=x^2的杆图'); subplot(2,2,4);

fill(x,y1,'r');%如果少了'r'则会出错 title('y1=x^2的填充图');

25

（3）>> x=-5:0.01:5; y=[];%起始设y为空向量 for x0=x

if x0<=0 %不能写成x0=<0

y=[y,(x0+sqrt(pi))/exp(2)]; %将x对应的函数值放到y中 else

y=[y,0.5*log(x0+sqrt(1+x0^2))]; end end plot(x,y)

26

（4）>> a=input('a='); b=input('b='); n=input('n='); t=-2*pi:0.01:2*pi; r=a*sin(b+n*t); polar(t,r) a=4 b=5 n=25

27

（5）>> x=linspace(-5,5,21); y=linspace(0,10,31);

[x,y]=meshgrid(x,y);%在[-5,5]*[0,10]的范围内生成网格坐标 z=cos(x).*cos(y).*exp(-sqrt(x.^2+y.^2)/4); subplot(2,1,1); surf(x,y,z); subplot(2,1,2); contour3(x,y,z,50);

（6） >>

ezsurf('cos(s)*cos(t)','cos(s)*sin(t)','sin(s)',[0,0.5*pi,0,1.5*pi]); %利用ezsurf隐函数

shading interp %进行插值着色处理

28

实验七 低层绘图操作

二（1）

h=figure('MenuBar','figure','color','r','WindowButtonDownFcn','disp(''Left Button Pressed'')')

h = 1

29

（2）>> x=-2:0.01:2; y=x.^2.*exp(2*x); h=line(x,y);

set(h,'color','r','linestyle',':','linewidth',2) text(1,exp(2),'y=x^2*exp(2*x)')

（3）>> t=0:0.00001:0.001; [t,x]=meshgrid(t);

30

v=10*exp(-0.01*x).*sin(2000*pi*t-0.2*x+pi); axes('view',[-37.5,30]); h=surface(t,x,v);

title('v=10*exp(-0.01*x).*sin(2000*pi*t-0.2*x+pi)'); xlabel(Ct'),ylabel('x'),zlabel('v') ??? Undefined function or variable 'Ct'.

（4）>> x=0:0.01:2*pi; y1=sin(x); y2=cos(x); y3=tan(x); y4=cot(x); subplot(2,2,1); plot(x,y1); subplot(2,2,2); plot(x,y2); subplot(2,2,3); plot(x,y3); subplot(2,2,4); plot(x,y4);

Warning: Divide by zero. > In cot at 13

31

（5）>> cylinder(5);

light('Position',[0,1,1]); material shiny

32

实验八 数据处理与多项式运算

（1）a >> A=rand(1,30000); b=mean(A) std(A,0,2) b =

0.5025 ans = 0.2899

B >> max(A) min(A) ans =

1.0000 ans =

1.5584e-005

C >> n=0; for i=1:30000 if A(i)>0.5 n=n+1; end end p=n/30000 p =

0.5042

（2） a >> A=45+51*rand(100,5); [Y,U]=max(A) [a,b]=min(A)

33

Y =

95.9123 95.6103 95.8956 95.6233 95.8496 U =

30 6 23 8 37 a =

45.1797 46.0073 45.5863 45.0149 45.1109 b =

51 46 29 21 17

B >> m=mean(A) s=std(A) m =

68.5490 66.6159 71.2299 71.4882 69.7301 s =

13.6861 14.3575 15.9196 14.4017 15.2025

C sum(A,2) [Y,U]=max(ans) [a,b]=min(ans) ans =

359.3754 373.0907 370.2539 342.8401 352.8497 425.8064

34

290.5902 391.0810 343.6866 341.2658 285.1690 383.5670 322.7536 290.7171 332.9745 351.3878 306.9964 374.7525 321.7565 296.2626 332.9247 351.9840 424.2966 327.2428 343.7518 374.2802 291.3005 421.8186 326.5739 405.9575 364.3370 374.2841 360.0279 289.0985 415.8491 297.5530 343.0118 394.9058 357.8329 348.0287 348.4106 344.1757 355.8473 339.9557 338.4695 342.0672 360.6158 330.5394 328.3173 369.8573

35

369.4287 336.8142 309.0138 415.6359 362.0495 337.1178 308.4595 357.7124 387.9822 359.9779 434.3310 399.2643 339.3533 330.8588 346.0471 323.5458 329.9120 328.8093 352.9168 403.3099 343.0933 346.9067 305.1942 309.4786 336.0630 320.6351 330.9037 344.8585 348.1910 297.4487 369.5013 289.9893 318.6062 288.7507 364.0922 371.7673 376.3963 356.6594 305.0560 342.5892 364.7775 362.2725 364.7526 357.5316

36

404.5707 396.4717 285.6594 350.2635 325.8875 289.9007 Y =

434.3310 U = 61 a =

285.1690 b = 11

C >> [zcj,xsxh]=sort(ans) zcj =

285.1690 285.6594 288.7507 289.0985 289.9007 289.9893 290.5902 290.7171 291.3005 296.2626 297.4487 297.5530 305.0560

37

305.1942 306.9964 308.4595 309.0138 309.4786 318.6062 320.6351 321.7565 322.7536 323.5458 325.8875 326.5739 327.2428 328.3173 328.8093 329.9120 330.5394 330.8588 330.9037 332.9247 332.9745 336.0630 336.8142 337.1178 338.4695 339.3533 339.9557 341.2658 342.0672 342.5892 342.8401 343.0118 343.0933 343.6866 343.7518 344.1757 344.8585 346.0471 346.9067 348.0287 348.1910 348.4106 350.2635 351.3878

38

351.9840 352.8497 352.9168 355.8473 356.6594 357.5316 357.7124 357.8329 359.3754 359.9779 360.0279 360.6158 362.0495 362.2725 364.0922 364.3370 364.7526 364.7775 369.4287 369.5013 369.8573 370.2539 371.7673 373.0907 374.2802 374.2841 374.7525 376.3963 383.5670 387.9822 391.0810 394.9058 396.4717 399.2643 403.3099 404.5707 405.9575 415.6359 415.8491 421.8186 424.2966 425.8064 434.3310

39

xsxh = 11 97 84 34 100 82 7 14 27 20 80 36 89 73 17 57 53 74 83 76 19 13 66 99 29 24 49 68 67 48 64 77 21 15 75 52 56 45 63 44 10

40

46 90 4 37 71 9 25 42 78 65 72 40 79 41 98 16 22 5 69 43 88 94 58 39 1 60 33 47 55 92 85 31 93 91 51 81 50 3 86 2 26 32 18 87

41

12 59 8 38 96 62 70 95 30 54 35 28 23 6 61

（3）>> h=6:2:18; x=6.5:2:17.5;

t1=[18,20,22,25,30,28,24]; t2=[15,19,24,28,34,32,30]; T1=spline(h,t1,x) T2=spline(h,t2,x) T1 =

18.5020 20.4986 22.5193 26.3775 30.2051 26.8178 T2 =

15.6553 20.3355 24.9089 29.6383 34.2568 30.9594

（4）>> x=1:0.1:101; y1=log10(x); p=polyfit(x,y1,5) y2=polyval(p,x);

plot(x,y1,':',x,y2,'-')

Warning: Polynomial is badly conditioned. Remove repeated data points or try centering and scaling as described in HELP POLYFIT. > In polyfit at 79 p =

42

0.0000 -0.0000 0.0001 -0.0040 0.1139 0.1641

（5）a >> x=1:0.1:101; y1=log10(x); p=polyfit(x,y1,5) y2=polyval(p,x);

plot(x,y1,':',x,y2,'-')

Warning: Polynomial is badly conditioned. Remove repeated data points or try centering and scaling as described in HELP POLYFIT. > In polyfit at 79 p =

0.0000 -0.0000 0.0001 -0.0040 0.1139 0.1641

>> p1=[1,2,4,0,5]; p2=[1,2]; p3=[1,2,3];

p=p1+[0,conv(p2,p3)] p =

1 3 8 7 11

43

B >> A=roots(p) A =

-1.3840 + 1.8317i -1.3840 - 1.8317i -0.1160 + 1.4400i -0.1160 - 1.4400i

C >> A=[-1,1.2,-1.4;0.75,2,3.5;0,5,2.5]; polyval(p,A) ans =

1.0e+003 *

0.0100 0.0382 0.0125 0.0223 0.0970 0.4122 0.0110 1.2460 0.1644

D >> polyvalm(p,A) ans =

1.0e+003 *

0.0076 -0.1281 -0.0775 0.1328 1.3900 1.1644 0.1824 1.7364 1.5198

实验九 数值微积分与方程数值求解

二（1）>> x=1; i=1;

f=inline('det([x x^2 x^3;1 2*x 3*x^2;0 2 6*x])'); while x<=3.01 g(i)=f(x); i=i+1;

x=x+0.01; %以0.01的步长增加，可再缩小步长提高精度end g;

t=1:0.01:3.01;

dx=diff(g)/0.01; %差分法近似求导

44

f1=dx(1) %x=1的数值倒数 f2=dx(101) %x=2的数值倒数 f3=dx(length(g)-1) %x=3的数值倒数 f1 =

6.0602 f2 =

24.1202 f3 =

54.1802

（2）>> f=inline('sqrt(cos(t.^2)+4*sin(2*t).^2+1)'); I1=quad(f,0,2*pi)

g=inline('log(1+x)./(1+x.^2)'); I2=quad(g,0,2*pi) I1 =

10.4285 I2 = 0.9997

（3）>> A=[6 5 -2 5;9 -1 4 -1;3 4 2 -2;3 -9 0 2]; b=[-4 13 1 11]'; x=A\\b y=inv(A)*b [L,U]=lu(A); z=U\\(L\\b) x =

0.6667

45

-1.0000 1.5000 -0.0000 y =

0.6667 -1.0000 1.5000 0.0000 z =

0.6667 -1.0000 1.5000 -0.0000

（4）function [x,y]=line_solution(A,b) [m,n]=size(A); y=[ ];

if norm(b)>0 %非齐次方程组 if rank(A)==rank([A,b]) if rank(A)==n disp('有唯一解x'); x=A\\b; else

disp('有无穷个解，特解x，基础解系y'); x=A\\b;

y=null(A,'r'); end else

disp('无解'); x=[ ]; end

else %齐次方程组 disp('有零解x'); x=zeros(n,1); if rank(A)

disp('有无穷个解，基础解系y');

46

y=null(A,'r'); end end

clc;clear; format rat

A=[2 7 3 1;3 5 2 2;9 4 1 7]; b=[6 4 2]';

[x,y]=line_solution(A,b) 运行结果

有无穷个解，特解x，基础解系y

Warning: Rank deficient, rank = 2, tol = 8.6112e-015. > In line_solution at 11 x =

-2/11 10/11 0 0 y =

1/11 -9/11 -5/11 1/11 1 0 0 1

（5）function g=f(x) g=3*x+sin(x)-exp(x);

clc;clear; fzero('f',1.5) ans =

1289/682

function F=fun(X) x=X(1); y=X(2); z=X(3);

F(1)=sin(x)+y^2+log(z)-7;

47

F(2)=3*x+2-z^3+1; F(3)=x+y+z-5;

X=fsolve('myfun',[1,1,1]',optimset('Display','off')) X =

909/1073 1735/728 1106/625

（6）function f=g(u) x=u(1); y=u(2);

f=2*x.^3+4*x.*y^3-10*x.*y+y.^2;

clc;clear; format long

f=inline('(x^3+cos(x)+x*log(x))/exp(x)'); [x,fmin1]=fminbnd(f,0,1)

[U,fmin2]=fminsearch('g',[0,0]) x =

0.522288340666172

fmin1 =

0.397363464998461 U =

1.001570135316681 0.833488282765738

fmin2 =

-3.324088491954234

（7）function xdot= sys( x,y) xdot=[y(2);(5*y(2)-y(1))/x];

clc;clear; x0=1.0e-9;xf=20;

48

[x,y]=ode45('sys',[x0,xf],[0 0]); [x,y]

x y’ y ans =

0.0000 0 0 0.5000 0 0 1.0000 0 0 1.5000 0 0 2.0000 0 0 2.5000 0 0 3.0000 0 0 3.5000 0 0 4.0000 0 0 4.5000 0 0 5.0000 0 0 5.5000 0 0 6.0000 0 0 6.5000 0 0 7.0000 0 0 7.5000 0 0 8.0000 0 0 8.5000 0 0 9.0000 0 0 9.5000 0 0 10.0000 0 0 10.5000 0 0 11.0000 0 0 11.5000 0 0 12.0000 0 0 12.5000 0 0 13.0000 0 0 13.5000 0 0 14.0000 0 0 14.5000 0 0 15.0000 0 0 15.5000 0 0 16.0000 0 0 16.5000 0 0 17.0000 0 0 17.5000 0 0 18.0000 0 0 18.5000 0 0 19.0000 0 0

49

19.5000 0 0 20.0000 0 0

（8）function xdot=sys(x,y)

xdot=[y(2)*y(3);-y(1)*y(3);-0.51*y(1)*y(2)];

clc;clear; t0=0;tf=8;

[x,y]=ode23('sys',[t0,tf],[0,1,1]) plot(x,y) x =

0 0.0001 0.0005 0.0025 0.0125 0.0625 0.1632 0.3033 0.4829 0.7162 0.9849 1.2610 1.5678 1.9550 2.3287 2.7024 3.0153 3.2921 3.4889 3.6452 3.7538 3.8624 3.9941 4.1645 4.3835 4.6537 4.9265 5.2245 5.5861 6.0302 6.3428 6.6555

50

本文来源：https://www.bwwdw.com/article/jzz8.html