计算机学院-课件-上机实验一 Matlab基础及科学计算

成绩:

Matlab基础及科学计算

姓名：黄晓阳

班号：HC001501

学号：2()1530()005

2345

4321

1、用Uatlab语句创建矩阵A,A=，如果给出“A(5,6)=你的学号”命令

6437

8946

将得到什么结果？

答:A方阵将扩充成5X6的矩阵,在第5行第6列处填充2015300005,多出来的其他未赋值处

填充0,原位置处为初始值.

源代码：

A=[2,3,4,5;4,3,2,1;6,4,3,7;8,9,4,6]

A(5,6)=2015300005

2、川magic命令生成一个六阶的魔方矩阵，并将其全部偶数行提取出来。

答:提取的结果为[8,30,4,32,28,36,6,2,34,26,22,12,10,14,18,24,20,16]

源代码：

Magic=reshape(magic⑹，[1,36]);

result=(];

fori=l:length(Magic)

ifrem(Magic(i),2)==0

result(length(result)+1)=Magic(i);

end

end

result

3、输入A为3x3的魔方阵，B为3x3的单位阵，由小矩阵组成3x6的大矩阵C和6x3

的大矩阵D,将D矩阵的最后1列构成向量E。

源代码：

A=magic(3)

B=ones(3,3)

C=zeros(3,6);

D=zeros(6,3);

fori=l:3

C(iz:)=[A(iz:),B(iz:)];

end

C

D=[A;B]

E=D(:,3)

123

4、输入矩阵A=456，使用全下标方式取出元素“3”，使用单下标方式取出元素

789

13

“8”，取出后2行子矩阵块,使用逻辑矩阵方式取出

79

答:详细操作见代码.

源代码：

A=[l,2,3;4,5,6;7,8,9]

three=A(1,3)

three_=A(7)

last_two_line=A(2:3,:)

A.*(rem(A,2)~=0).*(A~=5)

50

5、用数值方法求出5=,2'=1+2+4+…+2”+25。，试不采用循环的形式给出和式的数

r=0

值解。

答:5=2251799813685247

源代码：

symsk

symsum(2Ak,k,0,50)

6、试用循环结构找出“学号后三位”以下的所有质数。

答:[2,3,5]

源代码：

fori=l:005

if(isprime(i))

i

end

end

7、Fibonacci数列定义为4=4-+%一2/=3,4…为一个数值增长很快的数歹力

试编写程序寻找Fibonacci数列中第一个大于10000的元素

答:10946

源代码：

i=l;

whilefibonacci(i)<10000

i=i+l;

and

fibonacci(i)

8、使用while循环计算满足。(2灯-1)>学号后五位的最小，〃值及相应的积。

答:mmin=3乘积为：15

源代码：

clear

clc

i=l;

c=l;

whilec*(2*i-l)<005

c=c*(2*i-l);

i=i+l;

end

i

c*(2*1-1)

9、下面给出一个迭代模型

%|=1+券一14%

=0.3玉

写出求解该模型的M函数，如果取迭代初值为毛=0,%=0,那么请进行30000次迭代

求出一组x和y向量，然后在所有的'和以坐标处点亮一个点(注意不要连线)，最后绘制

出所需的图形。(注：绘制出的图形又称为，它将迭代出来的随机点吸引到一起，最后得出

貌似连贯的引力线图)

T9(30000)

functionT9(n)

x=zeros(size(1,n))；

y=zeros(size(1,n));

for1=1:n-l

x(i+l)=l+y(i)-1.4*x(i)A2;

y(i+l)=0.3*x(i);

end

plot(x,y,1.1)

title('Henon引力线图，)

10、绘制出下面的图形sin（；）,其中由于该函数在1=0附近震荡剧烈，因此,

尝试在在0附近进行调整步长。绘制并对比两种时间序列下的结果。第一个是采用0.02的

时间间隔，另外一个在/w[-0.2,0.2]附近将时间间隔设置为0.001。

*p，0.02step«0.001

源代码：

t=[-l:0.02:l];

subplot(1,2,1)

plot(t,oin(1./t))

title(*step=0.021)

t=-0.2:0.001:0.2;

subplot(1,2,2)

plot(t,sin(1./t))

title(*step=0.0011)

2*X2+1x>=\

11、绘出函数》二,0-Icxvl的图像。

-xx<=-1

源代码：

a=[-10:0.01:10];

plot(a(a>=l),2.*a(a>=l).A2+l,a(a<l&a>-l),0,a(a<=-l),-a(a<=-l).A3)

sin(86/3)

12、试绘制出p-的极坐标曲线，选择，变量的范围为(0,4%)和

学号尾号-cos?(3。/2)

(0,6幻，观察极坐标曲线的情况

源代码：

subplot(1,2,1)

theta=[0:0.01:4*pi];

rho=sin(8*theta/3)./(005-cos(3*theta/2).A2);

polar(theta,rho)

subplot(1,2,2)

theta=[0:0.01:6*pi];

rho=sin(8*theta/3)./(005-cos(3*theta/2).A2);

polar(theta,rho)

13、请分别绘制出p和sin。〉，）的三维图和等高线

10

源代码：

[x,y]=meshgrid(-10:0.1:10);

subplot(2,2,1)

meoh(x,y,x.*y)

subplot(2,2,2)

contour(x,y,x.*y,80)

subplot(2,2,3)

mesh(x,y,sin(x.*y))

subplot(2,2,4)

contour(x,y,sin(x.*y))

14、假设A为四阶希尔伯特矩阵，分别找出由矩阵A的第1、2、3、4条对角线上的元

素组成的列矢量。

-1

'1/2-

1/3~l/3~

答:diag_l=diag_2=1/4diag_3=diag_4=[l/4]

1/5_l/5_

1/6

1/7_

源代码：

A=hilb(4)

diag(A,0)

diag(A,1)

diag(A,2)

diag(A,3)

123

15、求矩阵A=3I2的转置矩阵、逆矩阵、矩阵的秩、矩阵的行列式值、矩阵的

_234

三次暴、矩阵的无穷范数，矩阵的特征值和特征向量。

-123~-132~909614「

答:A=312A'=2I3A3=8790132

_234_324132141207

-

2/31/31/3

inv(A)=8/3-2/1§7/3

7/31/3-5/3

7.290800

特征值矩阵为：（）-0.71680

00-0.5739

0.49720.17830.1051

对应的特征向量矩阵为：0.4691-0.8628-0.8520

0.72990.47310.5129

rank(A)=3|A|=3||A||-二9

源代码：

A=[l,2,3;3,l,2;2,3,4]

zhuanzhi=A'

ni=inv(A)

zhi=rank(A)

hanglieshi=det(A)

three_mi=AA3

wuqiongfanshu=norm(A,inf)

[tezhengxiangliang,tezhengzhi]=eig(A);

tezhengzhi

tezhengxiangliang

x+y+z=10

16、求方程组＜)，+z=0的解

2x-y-z=-4

x=2

答：，y=5

z=3

源代码：

[x,y,z]=solve(*x+y+z=10,x-y+z=O,2*x-y-z=-4','x,yzz*)

x}+x2-3X3-X4=1

17、求非齐次线性方程组卜%-3&+4七=4的通解

N+5X2-9X3-8X4=0

0"3/2--3/4

03/27/4

答：x=+G+C2

-8/1510

3/5_01

源代码：

A=[1,1,-3,-1;3,-1,-3,4;1,5,-9,-8];

b=[1;4;0];

B=[A,b];

teShuJie=A\b

tongjie=null(A,'r1)

18、求解下面的Lyapunov方程，并检验所得解的精度

0.31960.2360-0.2096

答：X:0.2624-0.18200.2844，将结果带入到所给矩阵运算式左端,可得右端为

-0.11780.5393-0.1797

1.000000000000003.000000000000002.00000000000000

3.000000000000004.00000000000000l.OOOOODOOOOOOOO，与上式右端相等，每

5.000000000000002.000000000000001.00000000000000

个元素的误差在10“4的量级.

源代码:

A=[l,2,3;4,5,6;7,8,0];

B=[2Z3,6;3,5Z2;4Z3Z1];

C=[l/3,2;3,4/l;5,2/L];

X=lyap(A,B,-C)

C_=A*X+X*B

C

A=8,.24k=0fl

19、计算zL，y——

.k\占Ax(A+l)

x

e'I

2

G:S,=A:(---)S,=l

X~X

源代码：

symskx

S2=symsum(1/(k*(k-1)),k,1,Inf)

A

SI=symsum(x(2*k)/factorial(k)zk,1,Inf)

20、对函数/=求函数的的一阶导和二阶导。

夕”

xcos(x)sin(x)xsin(x)

答：〈

xsin(x)cos(x)xcos(x)sin(x)

4

源代码：

symsx

yijiedao=diff(x*sin(x)/exp(2*x),1)

erjiedao=diff(x*sin(x)/exp(2*x),2)

21、人造地球卫星的轨迹可视为平面上的椭圆。我国第一颗人造地球卫星近地点距离地

球表面439km,远地点距离地球平面2384km,地球半径为6371km,求该卫星的轨迹长度。

人造地球卫星的轨迹可用椭圆的参数方程来表示，即

.r=«cosr

(0<t<2兀;a,b>0)

y=bsh\t

其中，a=6371+2384=8755(k⑼。=6371+439=6810(痴)。

由曲线积分可知，卫星估计长度L为积分L=4亢°2si/r+b28s2dt的值。

答:该卫星轨迹长度为4.9095e+04(km).

源代码：

%T21(2384+6371,439+6371,10000)

functionT21(a,b,n)

len=0;

step=pi/2/n;

fortheta=0:step:pi/2

len=len+Len(a*cos(theta),a*cos(theta+stap),...

b*sin(theta),b*sin(theta+step));

and

perimeter=len*4

end

functiony=Len(xlzx2,yl,y2)

y=((xl-x2)A2+(yl-y2)A2)A0.5;

end

22、求二重积分J'JJ(ysinx+xsi”)dMv，]了,—+7)，+6"dt

答:S]=-98.6960S2=14.5

源代码：

symsxy

f1=@(xzy)y.*sin(x)-x.*sin(y);

f2=@(x,y)xJ2+7.*yJ2+6;

integral2(fl,3*pi,4Tpi,pi,2*pi)

integral2(£2,0,1,@(x)x,@(x)x+1)

23、求解微分方程的通解和初始条件为M0=0，W0=l时的特解，并采用ode45()的数

值求解方法解该方程组，绘图并比较两种求解方法。

公

-

小

电

力

ode<45)求杆x(t)用像

源代码：

functionT23()

[x,y]=dsolve('Dx=y+x','Dy=2*x1,'t*);

disp(x)

disp(y)

[x,y]=dsolve('Dx=y+x*,'Dy=2*x1,'x(0)=0*,'y(0)=1*,'t*);

disp(x)

disp(y)

subplot(2,2,1)

ezplot(x)

title(*dsolve求解x(二)图像，)

gridon

axis([0,4,0,1000])

subplot(2,2,2)

ezplot(y)

title('dsolve求解y(二)图像，)

gridon

axis([0z4,0z1000])

g=@(t,y)Wf(t,y);

[t,x]=ode45(g,[0,4.2],[0,1]);

subplot(2,2,3)

plot(t,x(:,1));

title(•ode(45)求解x(t)图像,)

gridon

axis([0z4,0z1000])

subplot(2,2,4)

plot(t,x(:,2))

title('ode(45)求解y(t)图像，)

gridon

axis([0z4z0,1000])

end

functiondy=Wf(t,y)

dy=zeros(2,1);

dy(l)=y(2)+y(l);

dy(2)=2*y(l);

end

24、求解下面的最优化问题

,，牒“(xTx+0

兑％=1

答:《

x2=0

源代码：

functionT24()

fmincon(@(x)x(1)A2-2*x(1)+x(2),[50,300]，匚,[],[],[],[],[],@nonl

con)

end

function[C,Ceq]=nonlcon(x)

C(l)=4*x(l)A2+x(2)A2-4;

C(2)=-x(l);

C(3)=-x(2);

Ceq=[];

end

25、求解/(1)=10012-4『+(1-西)2在1121]上的最小值

答:当未知量都取1时，有最小值f(X)min=0

源代码：

functionT25()

fmincon(@(x)100*(x(2)-x(l)A2)A2+(l-x(l))A2,[0,0],[],[],[],[],[J,[：/

un)

end

function[c,ceq]=nun(x)

c(l)=-x(l)-1.2;

c(2)=x(l)-l;

c(3)=-x(2)-1.2;

c(4)=x(2)-1;

ceq=[];

end

min(X)+3/-5x3)

的线性规划

26、计算4

3%+鼻W4

xpx2,x3>0

尤i=0

答:当未知量取值为,占=0时，有最小值-20

£=4

源代码：

functionT26()

fmincon((a(x)x(l)+3*x(2)-5*x(3)/[0,0,0],[],：],[],[],[],[],@nun)

end

function[c,ceq]=nun(x)

ceq=[];

c(1)=x(1)-x(2)+x(3)-8;

c(2)=3*x(l)+x(3)-4;

c(3)=-x(1);

c(4)=-x⑵;

c(5)=-x(3);

end

27、已知x=0:0.2:2,函数y=(x?-3x+学号尾号)，"sinx。对外=0:0.03:2,分别采

用最近邻插值、后临插值、线性插值、三次样条插值方法进行插值，并画出响应的图形。

近邻插僮E邻插俶

源代码：

clea