中南大学数学建模报告

来源：华拓科技网

数学实验与数学建模实验报告

学院：信息科学与工程学院
专业班级：
姓名：
学号：

完成时间： 2014年1月1日

承诺书

本人承诺所呈交的数学实验与数学建模作业都是本人通过学习自行进行编程完成，所有结果都通过上机验证，无转载或抄袭他人，也未经他人转载或抄袭。若承诺不实，本人愿意承担一切责任。

承诺人：边浩然

2014年1月6日
注意事项如下：
1、上机时间：第11周到第18周星期六晚上：7：00-8：30；
2、上机地点：新校区数学与统计学院二楼数学实验室；3、交作业时间：第19周的星期五（2014年1月6日）交到新校区数学与统计学院二楼数学实验室数学实验室办公室；4、报告所有结果交电子打印稿到新校区数学与统计学院二楼数学实验室办公室，并将电子文档发送到邮箱：（word文档命名：姓名＋学号＋数学实验作业）

数学实验学习体会
（每个人必须要写1500字以上，占总成绩的20%）

大二上学期的时候，我选修了数学建模这门课，一开始以为数学建模这门课主要是对数学方面的学习。所以上第一节课的时候很迷茫，因为第一节课都讲的是MATLAB这个软件的运用。与去年的高等数学，概率论与数理统计和线性代数相比，数学建模实验除了有对数学基础的要求外，还要我们有一定的编程能力，尽管不像C语言那样有严格的语法，但也有不算简单的编程规则。事实证明，选这门课对我学习掌握MATLAB是十分有帮助的！

在上理论课的时候对MATLAB这个软件的第一认识是，这是一个类似于C语言编程的编辑软件，使用软件的主要手段就是编写代码，通过编写代码来告诉软件需要制作什么样的模型，所以第一件事就是记代码的编写规则。由于大一上学期有C语言的基础，对于代码的使用与学习都有一定的优势，在学MATLAB中的各种表达式都有很比较容易。但由于MATLAB 有一部分是很细小的差别，所以在记忆的中有一部分的编写规则与C语言相差还是很大的，
时候总是容易记混。但好在没有特别大的问题，不会对整个学习过程有大影响。

纯理论的知识与经过实践后的知识的无法匹比的，像MATLAB这样注重实践的课更是这样，所以这次我们上完十六个学时的理论课后，还有十六个学时的实验课。为的就是让我们更加充分的掌握这么软件。为的就是让我们更加充分的掌握这么软件。在刚刚开始到实验室去做实验的时候，我明显有一些不适应，一是因为着实被这么多的题目吓了一跳，匆匆浏览了一下，惊恐地发现大部分貌似都不太会，当时就感觉前途一片迷茫，二是因为对Matlab软件不熟悉，另外就是对于之前一年所学的数学知识的记忆已经模糊了。但我仔细一想，后面还有两个月的时间，只要慢慢来，有条不紊，一定可以做完这些题的。于是我开始理论复习与实验相结合地完成作业，从图书馆借了一本有关MATLAB的书籍，翻阅相关的代码还有例题。

数经过大概一个月左右的实验，我发现一开始我对数学建模这个课程的理解是错误的。

学建模这个课程不是仅仅教授我们怎么样使用MATLAB这个软件，更注重的是数学思想上的培养，因为在实验中会遇到很多数学函数，有图型的函数，有导数的函数等等。在做MATLAB我基本上是按实验顺序和题目顺序来做，如果偶尔碰到不会做的我会首先查实验的过程中，阅有关的书籍，看看书上有没有类似的解决办法。如果有，那么我会仿照书上的方法来解决这个题目；如果没有，那么我会再查看有关课件或者询问同学来完成此题目的解答。但其实很多题即使是这样也还是不能解决的，因为数学知识有限或者说算法思想不正确导致当时不

能理解错误在哪，不过后来做了很多其他的题目中，有时可以渐渐领悟一些东西，让我恍然大悟，突然想起以前一道没做出来的题该怎么做。这种在实验里领悟到的经验、体会都是无可替代的，它们才是在实际运用中占据主要地位的知识。

在做完试验后我对数学建模认识发生的一定的改变，也更加深刻的了解了这门课程，数理论科学方面有着非常学建模的是一门基于数学理论基础的实践性科学，数学建模在工科、
多的应用，与我们电学紧密相关的比如“小波分析”，另外还有接下来我们下个学期要学的“系，“系统仿真技术”就是关于matlab在电气工程专业上的应用。既然数学建模统仿真技术”
是一门实践性科学，我们就要不断地去练习，去实践，比较不同算法的相异之处，若想要完成这些，我们就要有着扎实的理论基础，而数学建模的理论实践基础就是建立在我们之前学过的微积分、线性代数等数学课程上，但是就我们现在所学的数学知识来看不够，在学习中如果想要真正学好数学建模，我我就发现数学建模中运用到了很多我们没学过的数学知识，
们还要继续更深入的学习数学，数学是自然科学的基础，只有掌握好数学，才能够在自然科学的道路上走的更远。

这学期的数学建模课程，让我收获良多，不但认识到MATLAB这个软件的方便与强大，也复习了很多数学方面的知识，同时我发现了自己在数学知识上的欠缺，并且也改正了不少。经过这样自己动手，通过理论与实践的结合，这样所掌握的知识与能力，将会比在课堂上，通过老师讲授所掌握的要深刻得多，经过实践，我们获得的是具象，而且有深刻的印象。这也是此次实验带给我的除了实践能力之外的收获。

实验一图形的画法

1.做出下列函数的图像：

（1）

(

)

sin(

2 )

，

（分别用plot、fplot）

x=-2:0.1:2;y=(x.^2).*sin(x.^2-x-2)；
fplot(x,y)；

x=-2:0.5:2;y=x.*x.*sin(x.^2-x-2);
plot(x,y)

（2）

/ 25

?（用参数方程）

t=-pi:pi/100:pi;x=3*cos(t);y=5*sin(t);plot(x,y);

(3) 在同一图形窗口中，画出四幅不同图形（用subplot命令）：

y 1

cos( )

，

sin(

/ 2)

，

y 3

cos(

)

，

sin( )

（

[

0 ,

2?]

）

x=linspace(0,2*pi,60);
y=cos(x);z=sin(x-0.5*pi);t=(x.^2).*cos(x-pi);ct=exp(sin(x));subplot(2,2,1);
plot(x,y);
title('cos(x)');axis([0,2*pi,-1,1]);
subplot(2,2,2);
plot(x,z);
title('sin(x-0.5*pi)');axis([0,2*pi,-1,1]);
subplot(2,2,3);
plot(x,t);
title('(x.^2).*cos(x-pi)');axis([0,2*pi,-40,40]);
subplot(2,2,4);
plot(x,ct);
title('exp[sin(x)]');axis([0,2*pi,-1,1];

2 作出极坐标方程为

2 ( 1

cos

)

的曲线的图形.

theta=0:0.01:2*pi;rho=2*(1-cos(theta));polar(theta,rho)；

3 作出极坐标方程为

r ?

e t

的对数螺线的图形.

t=0:0.01:6*pi;
rho=exp(t/10);
polar(t,rho,'k');
title('对数螺线的图形');

4 绘制螺旋线

???

cos

在区间［

，

］上的图形.在上实验中，显示坐标轴名称

sin

t=0:0.1:4*pi;x=4*sin(t);
y=4*cos(t);z=t;plot3(x,y,z,'r');
gridon;xlabel('x_axis');
ylabel('y_axis');zlabel('z_axis');
title('空间曲线');

5 作出函数

xye

的图形.

x=-2:0.1:2;y=x;
[X,Y]=meshgrid(x,y);
Z=-X.*Y.*exp(-X.^2-Y.^2);

mesh(X,Y,Z);

6 作出椭球面	x	2	?	y	2	?	z	2	?	1	的图形.
	4			9			1

(该曲面的参数方程为

sin

cos

sin

cos

(

?, 0

2?).)

u=0:0.1:pi;v=0:0.1:2*pi;[u,v]=meshgrid(u,v);x=2*sin(u).*cos(v);y=3*sin(u).*sin(v);z=cos(u);
mesh(x,y,z)

7 作双叶双曲面	x	2		?	y	2		?	z	2		?	?1	的图形.
	1 . 5		2		1 . 4		2		1 . 3		2

(曲面的参数方程是

1 . 5

cot

cos

1 . 4

cot

sin

1 . 3

csc

0 ,

?时对应双叶双曲面的另

其中参数

时对应双叶双曲面的一叶, 参数

一叶.)
ezmesh('1.5*cot(u).*cos(v)','1.4*cot(u).*sin(v)','1.3*csc(u)',[pi/1000000,pi/2,-pi,pi]);hold on;
ezmesh('1.5*cot(u).*cos(v)','1.4*cot(u).*sin(v)','1.3*csc(u)',[-pi/2,pi/1000000,-pi,pi]);axis([-15,15,-15,15,-15,15])

8作出圆环

( 8

cos

)

cos

( 8

cos

)

sin

3?/

2 ,

)

的图形.

ezmesh('(8+3*cos(u))*cos(v)','(8+3*cos(u))*sin(v)','7*sin(u)',[0,2*pi,0,2*pi]);axis([-20,20,-20,20,-20,20])
title('圆环的图形')

9 作出球面

和柱面

(

1 )

相交的图形.

symsu v;
u=0:0.2:2*pi;
v=0:0.2:2*pi;
[u,v]=meshgrid(u,v);
x1=2*sin(u).*cos(v);
y1=2*sin(u).*sin(v);
z1=2*cos(u);
x2=sin(u)+1;
y2=cos(u)/sqrt(2);
z2=sin(v);

mesh(x1,y1,z1);
holdon;
surf(x2,y2,z2);
holdoff;

10 作出锥面

和柱面

(

1 )

相交的图形.

symsu v;
u=0:0.2:2*pi;
v=0:0.2:2*pi;
[u,v]=meshgrid(u,v);
x1=sin(u).*sin(v);
y1=sin(u).*cos(v);
z1=sin(u);
x2=sin(u)+1;
y2=cos(u);
z2=sin(v);
mesh(x1,y1,z1);
holdon;
surf(x2,y2,z2);
holdoff;

11用动画演示由曲线y?sinz,z?[0 ,?]绕z轴旋转产生旋转曲面的过程.（该曲线绕z轴旋转所得旋转曲面的

方程为x2?y2?sin2z,其参数方程为

x?sinzcosu,y?sinzsinu,z?z, (z?[0 ,?],u?[0 , 2?])）

>> m=moviein(20);

fori = 1:20

u=0:0.1:pi/5*(i+0.2);

v=0:0.1:pi;

[u,v]=meshgrid(u,v);

x=sin(v).*cos(u);

y=sin(v).*sin(u);

z=v;

mesh(x,y,z)

m(:,i)=getframe;

end

movie(m,2)

12. 画出变上限函数	?	x	t	sin	t	2	dt	及其导函数的图形.
		0

symst;

x=-2*pi:0.1:2*pi;

int(t*(sin(t))^2,0,x);
y1=(x.*(sin(x)).^2);
y2=(-1/2*x.*cos(x).*sin(x)+1/4*x.^2+1/4*sin(x).^2);plot(x,y1,x,y2)

13.迪卡尔曲线

3 at

(

3 axy

0 )

t=-pi:pi/50:pi;
x=(3*t)./(1+t.*t);
y=(3*t.*t)./(1+t.*t);
plot(x,y)

14.蔓叶线

(

)

>>t=-pi:pi/50:pi;
x=(t.*t)./(1+t.*t);
y=(t.*t.*t)./(1+t.*t);
plot(x,y)

15.摆线

( t

sin

( 1

cos

)

t=-pi:pi/50:pi;
x=t-sin(t);
y=1-cos(y);
plot(x,y)

																2			2			2
16.内摆线（星形线）	x	?	a	cos	3	t	,	y	?	a	sin	3	t	(	x	3	?	y	3	?	a	3	)

t=-pi:pi/50:pi;
x=(cos(t).^3);
y=(sin(t).^3);
plot(x,y)

17.圆的渐伸线（渐开线)

(cos

sin

(sin

cos

)

t=-pi:pi/50:pi;
x=cos(t)+t.*sin(t);
y=sin(t)-t.*cos(t);
plot(x,y)

18.空间螺线

cos

sin

t=-pi:pi/50:pi;
x=cos(t);
y=sin(t);
z=t;
plot3(x,y,z)
gridon;xlabel('x_axis');
ylabel('y_axis');zlabel('z_axis');
title('空间曲线');

19.阿基米德线

。

t=-pi:pi/50:pi;
x=t;
polar(t,x)

20.对数螺线

r ?

。

t=-pi:pi/50:pi;
x=exp(t);
polar(t,x)

21.双纽线

cos

2?((

)

(

))

t=-pi:pi/50:pi;
r=sqrt(cos(2*t));
polar(t,r)

22.双纽线

sin

2?((

)

t=-pi:pi/50:pi;
>>r=sqrt(sin(2*t));
>>polar(t,r)

23.四叶玫瑰线

sin

t=-pi:pi/50:pi;
r=sin(2*t);
polar(t,r)

24.玫瑰线

sin

t=-pi:pi/50:pi;
r=sin(3*t);
polar(t,r)

25.三叶玫瑰线

cos

t=-pi:pi/50:pi;
r=cos(3*t);
polar(t,r)

26.作出以参数方程表示的空间曲线

0 . 2

cos

?
t

?
e

0 . 2

sin

[

0 ,

20 ]

t=0:0.1:20;
x=exp(0.2*t).*cos((pi/2).*t);
y=(pi/2)*exp(0.2*t).*sin(t);
z=t;

plot3(x,y,z);
gridon;xlabel('x_axis');
ylabel('y_axis');zlabel('z_axis');

27.以绘制极坐标系下曲线

cos( b

n?)

，并讨论参数

b ,

的影响。

t=-pi:pi/50:pi;
subplot(2,2,1);
polar(t,cos(t));
title('a=b=n=1');
subplot(2,2,2);
polar(t,10*cos(t));
title(a=10,b=n=1);
title('a=10,b=n=1');
subplot(2,2,3);
polar(t,cos(10+t));
title('b=10,a=c=1');
subplot(2,2,4);
polar(t,cos(10*t));
title('n=10,a=b=1')

28. (曲线族绘制) 三次抛物线的方程为

，试探讨参数a 和c 对其图形的影响。

x=-2:0.1:2;
subplot(2,2,1);
plot(x,x.^3+x)
title('a=c=1');
subplot(2,2,2);
plot(x,5*x.^3+x)
title('a=5,c=1');
subplot(2,2,3);
plot(x,x.^3+5*x)
title('a=1,c=5');
subplot(2,2,4);
plot(x,5*x.^3+5*x)
title('a=c=5');

sin

30.画出空间曲线

在

范围内的图形，并画出相应的等高线。

[a,b]=meshgrid(-30:1:30);
c=sqrt(a.^2+b.^2)+eps;
d=1+sqrt(a.^2+b.^2)+eps;
z=sin(c)./d;
mesh(a,b,z)
contour(a,b,z)

31.根据给定的参数方程，绘制下列曲面的图形。

椭球面

cos

sin

cos

sin

u=-pi:pi/50:pi;

v=-pi:pi/50:pi;
[u,v]=meshgrid(u,v);
x=3*cos(u).*sin(v);
y=2*cos(u).*cos(v);
z=sin(u);
mesh(x,y,z)

椭圆抛物面

3 u

sin

2 u

cos

4 u

u=-pi:pi/50:pi;
v=-pi:pi/50:pi;
[u,v]=meshgrid(u,v);
x=3*u.*sin(v);
y=2*u.*cos(v);
z=4*u.*u;
mesh(x,y,z)

单叶双曲面

sec

sin

sec

cos

tan

u=-pi:pi/50:pi;

v=-pi:pi/50:pi;
[u,v]=meshgrid(u,v);
x=3*sec(u).*sin(v);
y=2*sec(u).*cos(v);
z=4*tan(u);
mesh(x,y,z)

双曲抛物面

u=-pi:pi/50:pi;
v=-pi:pi/50:pi;
[u,v]=meshgrid(u,v);
x=u;
y=v;
z=(u.*u-v.*v)./3;
mesh(x,y,z)

旋转面

sin

cos

u=0:0.1:5;
v=0:0.1:5;
[u,v]=meshgrid(u,v);
x=log(u).*sin(v);
y=log(u).*cos(v);
z=u;
mesh(x,y,z)

圆锥面

sin

cos

u=-pi:pi/50:pi;
v=-pi:pi/50:pi;
[u,v]=meshgrid(u,v);
x=u.*sin(v);
y=u.*cos(v);
z=u;
mesh(x,y,z)

环面

( 3

0 . 4

cos

)

cos

( 3

0 . 4

cos

)

sin

0 . 4

sin

>>u=-pi:pi/50:pi;

>>v=-pi:pi/50:pi;
>>[u,v]=meshgrid(u,v);
>>x=(3+0.4*cos(u)).*cos(v);
>>y=(3+0.4*cos(u)).*sin(v);
>>z=0.4*sin(v);
>>mesh(x,y,z)

正螺面

sin

cos

4 v

>>u=-pi:pi/50:pi;
>>v=-pi:pi/50:pi;
>>[u,v]=meshgrid(u,v);
>>x=u.*sin(v);
>>y=u.*cos(v);
>>z=4*v;
>>mesh(x,y,z)

实验二一元函数微分学

1.分别画出坐标为

( i

4 i

( i

1 ,

2 ,

, 10 )

的散点图, 并画出折线图.

>>for i=1:10
plot(i,i.^2,'.');
hold on
plot(i.^2,4*i.^2+i.^3,'.');
end

>>x=1:10;
>>y=x.^2;
>>plot(x,y);
>>plot(x.^2,4*x.^2+x.^3);
>>axis([0,105,0,1450])

2.画出前25个素数的散点图.

i=1;
j=1;
prinum=zeros(25,1);
while(i<=25)
if(isprime(j))
prinum(i)=j;
i=i+1;
end
j=j+1;
end
plot(prinum,'o')
holdon;

3. 设数列

{ n

}

与

{ n

}

由下式确定:

1 ,

(

)

观察

{ n

}

与

{ n

}

的极限是否存在.

f[x_,y_]:=Sqrt[x*y];
g[x_,y_]:=(x+y)/2;
xn=1;yn=2;
For[n=1,n<=10,n++,xN=xn;yN=yn;xn=N[f[xN,yN]];yn=N[g[xN,yN]];Print["x",n,"=",xn," y",n,"=",yn]]

4.讨论极限

lim

n??

cos

>>syms x n
>>limit((cos(x)).^n,n,inf)

ans=

piecewise([cos(x)= 1, 1], [1 < cos(x) and 1 < abs(cos(x)), Inf], [not 1 <=cos(x) and 1 < abs(cos(x)), limit(cos(x)^n, n = Inf)],[abs(cos(x)) in (0, 1) and cos(x) <> 1, 0], [(1 < cos(x) orcos(x) in (0, 1)) and abs(cos(x)) = 1, NaN])
5.在MATLAB中求下列极限（写出MATLAB命令和运行结果）

(1)

lim

n??

(

)

（2）

lim

x??

( 1

)

（3）

lim x?0

sin

（4）

lim x??

（5）

lim

n??

(6)

lim

n???

(

?1 )

5 n

3 n

(7)

lim x?0

sin

cos

(8)

lim

x??0

cot

（9）

lim
x??0

sin

（12）

lim?

x?0?

sin

(10)

lim
x??0

（11）

lim

x??

(sin

cos

)

cos

（1）symsn
limit(sqrt(n+sqrt(n))-sqrt(n),n,inf)
ans =1/2
（2）symsx
limit((1-2/x).^(3*x),x,inf)
ans=1/exp(6)
（3）symsx
limit(sin(x)./(x.^3+3*x),x,0)
ans =1/3
(4)syms x
limit((3*x.^3-4*x.^2+2)/(7*x.^3+4),x,inf)
ans =3/7
(5)symsn

limit((2*n.^3+1)/(5*n.^3+1),n,inf)
ans=2/5
(6)syms n
limit(((-1).^n+4.^n)/(3.^(n+1)+4.^(n+1)),n,inf)ans=1/4
(7)symsx;
limit((sin(x)-x*cos(x))/x^2/sin(x),x,0)
ans=1/3
(8)symsx;
limit(log(cot(x))/log(x),x,0,'right')
ans=-1
(9)symsx;
limit(x^2*log(x),x,0,'right')
ans=0
(10)symsx;
limit(x^x,x,0,'right')
ans=1
(11))syms x
limit((sin(1/x)+cos(1/x)).^x,x,inf)
ans=exp(1)
(12)symsx;
limit((sin(x)/x)^(1/(1-cos(x))),x,0)
ans=exp(-1/3)=0.7165
6.讨论下列函数在指定点的连续性：

函数

(

)

???

x?1

在

处的连续性；

symsx y;
y=(x.^2-5*x-6)./(x+1);
a=limit(y,x,-1,'left');
b=-7;
ifa==b;
disp('该函数在该点有极限！函数极限为：'); a
else
disp('该函数在该点无极限！');
end
该函数在该点有极限！函数极限为：a=-7

函数

(

)

??1

在

处的连续性；

sin

symsx y;
y1=1-x.^2-5*x;

y2=sin(x)./x;
a=limit(y1,x,0,'left');
b=limit(y2,x,0,'right');
c=1;
ifa==b& b==c;
disp('该函数在该点有极限！函数极限为：'); a
else
disp('该函数在该点无极限！');
end
该函数在该点有极限！函数极限为：a=1

7.根据要求在MATLAB中求下列函数的导数

，求

f x ( )

arcsin

?
?
?，求

(1)

(2)

f ?

（3）设

，求dy(4)

ln(1

)

，求

dx 2

x?1

(1)>>syms a x y;y=a^a+a^x+x^a+x^(a*x);diff=diff(y,x)
diff=a*x^(a - 1) + a^x*log(a) + a*x*x^(a*x - 1) + a*x^(a*x)*log(x) (2)>>syms x;y=diff(asin((1-x^2)/(1+x^2)),x),x=1;eval(y)
y=(-2*x/(1+x^2)-2*(1-x^2)/(1+x^2)^2*x)/(1-(1-x^2)^2/(1+x^2)^2)^(1/2) ans=-1
(3)>>syms x
>>dy=diff(log(x+sqrt(a^2+x^2)))
dy=
(1+1/(a^2+x^2)^(1/2)*x)/(x+(a^2+x^2)^(1/2))
(4>>syms x
>>y=x^2*log(1+x);
>>diff(y,2),x=1;eval(ans)
ans=
2*log(1+x)+4*x/(1+x)-x^2/(1+x)^2
ans =
2048/653

8.函数		f	(	x	)	?	1	/	x	4	在区间[1,2]上满足拉格朗日中值定理的条件, 因此存在					??	( 1 ,	2 )	使
f ?)	?	f	(	f	(	2 )	?		f	( 1 ))		/(	2	?	1 ).	??	( 1 ,	2 )	使

可以验证这个结论的正确性.

>>clear;
clc;
y= sym('4*x^3-5*x^2+x-2'); %%%% y
y0= subs(y,'x',0); %%%%% y(0)
y1= subs(y,'x',1); %%%%% y(1)

d_y= diff(y); %%%%% y'
kesi= solve(d_y*(1-0)-(y1-y0)); %%%%% 解y'(ξ)*(b-a)=f(b)-f(a)%%%%% 接下来需要判断ξ是否是[0,1]之间 %%%%%%%%%%%
kesi= double(kesi)
forii = 1:length(kesi)
if kesi(ii)>0&kesi(ii)<1
disp('ξ∈[0,1],Lagrange mean value theorem is proved!'); break;
end
end
kesi=

0.11620.7171ξ∈[0,1],Lagrange mean value theorem is proved!

9.验证拉格朗日定理对函数

在区间[0,1]上的正确性.

y= sym('4*x^3-5*x^2+x-2');
y0= subs(y,'x',0);
y1= subs(y,'x',1);
d_y= diff(y);
kesi= solve(d_y*(1-0)-(y1-y0));
kesi= double(kesi)
forii = 1:length(kesi)
if kesi(ii)>0&kesi(ii)<1
disp('ξ∈[0,1],Lagrange mean value theorem is proved!'); break;
end
end

10.证明:对函数

应用拉格朗日中值定理时, 所求得的点

?总是位于区间

[

]

的正中间.

symsx p q r a b;y=p*x^2+q*x+r;diff(y) ans=2*p*x+q
f=(2*p*x+q)*(a-b)-((p*a^2+q*a+r)-(p*b^2+q*b+r));c=solve(f) c=1/2*a+1/2*b

11.

ln( 1

)

，求

dx 2

x?1

>>syms x
>>y=x.^2*log(1+x);
>>dy=diff(y,2);
>>x=1;
>>eval(dy)
ans=3.1363

12.设	??	x	?	a	( t	?	sin	t	)	，求	dy

	?	y	?	a	( 1	?	cos	t	)		dx

>>syms t
>>syms a
>>x=a*(t-sin(t));
>>y=a*(1-cos(t));
>>dy_dx=diff(y,t)/diff(x,t)
dy_dx =-sin(t)/(cos(t) - 1)

13．已知多项式

f x ( )

?，

g x ( )

，求：

（1）

)

的根； (2)

)

在闭区间[-1,2]上的最小值；

(

)

（3）

(

)

(

)

，

(

)

(

)

和

(

)

；（4）

)

的导数。

(1)
>>syms x
>>f=6*x.^5+2*x.^3-5*x.^2+1;
>>solve(f,x)

ans=

-0.4018369413802212127835619534773
0.091486951782804328430730634259985171062112272218481079452136829*i +0.605736693623540593181038577821345292393238255313202685865465

0.605736693623540593181038577821346717573441928475188877857986 -0.09148695178280432843073063426*i
-
0.96851301346816273468061031925246*i- 0.40557110297925174766753712711913

0.96851301346816273468061031925245840435326955993523868688442211*i -0.40557110297925174766753712711913
(2)
（3）>>syms x
>>f=6*x.^5+2*x.^3-5*x.^2+1;
>>g=(1/6)*x.^4+2*x.^3-5*x.^2+1;
>>h=f+g
h =6*x^5 + x^4/6 + 4*x^3 - 10*x^2 + 2
>>h=f*g
h =(x^4/6 + 2*x^3 - 5*x^2 + 1)*(6*x^5 + 2*x^3 - 5*x^2 + 1)
>>h=f/g
h =(6*x^5 + 2*x^3 - 5*x^2 + 1)/(x^4/6 + 2*x^3 - 5*x^2 + 1)
(4)>> diff(f)

ans=30*x^4 + 6*x^2 - 10*x
14.已知函数

f	(	x	)	?	1	x	6	?	2	x	5	?	25	x	4	?	60			3	?	15		x	2	?	180		x	?	25 ,
					2								2							f							180
在区间					[?6 , 6 ]				上画出函数							f	(	x	),	f	?(	x	),	f	??(	x	)	的图形, 并找出所有的驻点和拐点.

>>syms x
>>f=(1/2)*x.^6-2*x.^5-(25/2)*x.^4+60*x.^3-150*x.^2-180*x-25; >>x=-6:0.1:6;
>>subplot(2,2,1);
>>plot(x,(1/2)*x.^6-2*x.^5-(25/2)*x.^4+60*x.^3-150*x.^2-180*x-25) >>title('f(x)');
>>diff(f,2)
ans=15*x^4 - 40*x^3 - 150*x^2 + 360*x - 300
>>subplot(2,2,2);
>>plot(x,15*x.^4 - 40*x.^3 - 150*x.^2 + 360*x - 300);
>>title('df(x)');
>>subplot(2,2,3);
>>diff(f,3)
ans=60*x^3 - 120*x^2 - 300*x + 360
>>plot(x,60*x.^3 - 120*x.^2 - 300*x + 360);
>>title('(df(x))^2')

15.求函数

sin

)

的位于区间

( ?)

内的极值的近似值.

functiony=f(x)
y=2*sin(2*x)*sin(2*x)+5/2*x*cos(x/2)*cos(x/2);ezplot(y,[0,pi]);
grid;
x=fminbnd('f1(x)',0.5,2.5)
f1(x)
x=fminbnd('-f1(x)',0,pi)
f1(x)

x=fminbnd('-f1(x)',1.5,pi)
f1(x)
极小值点x= 1.6239 ans = 1.9446
极大值点x= 0.82 ans = 3.7323
极大值点x= 2.2449 ans = 2.9571

实验三一元函数积分学

一元函数积分学
1．用MATLAB计算下列不定积分。

（1）

（2）

sin cos

xdx

（1）
>>syms x
>>f=sqrt(x.^2+1)/x.^2;
>>int(f)
ans=asinh(x) - (x^2 + 1)^(1/2)/x
（2）
>>syms x a;
>>f=a.^x*sin(x).*(cos(x)).^2;
>>int(f)
ans=(a^x*(log(a)^3*cos(x)^2*sin(x)-log(a)^2*cos(x)^3+2*log(a)^2*cos(x)*sin(x)^2 +3*log(a)*cos(x)^2*sin(x)+2*log(a)*sin(x)^3-3*cos(x)^3))/(log(a)^4+10*log(a)^2 + 9)
2．用MATLAB求解下列各积分。

（1）

cos

xdx

（2）

sin 2 tdt

。

（3）设

f x ( )

???

，求

2?0

f x dx

(1)>> syms x;
>>v=int(exp(2*x).*cos(x),0,2*pi)

v=(2*exp(4*pi))/5 - 2/5
>>vpa(v)
ans=114700.1252546613198676990947
(2)>>syms t;
v=int(exp(-t)*sin(2*t),0,inf)
v = 2/5
>>vpa(v)
ans=0.4
(3)>> syms x
>>v=int(x.^2,0,1)+int(x,1,2)
v =11/6

4．求由曲线

(

绕x 轴旋转所产生的旋转体的体积。

symsx;y=pi*(5+sqrt(16-x.^2)).^2;int(y,x,-4,4); ans=856/3*pi+80*pi^2
5．求下列曲线与所围成图形的面积：

（1）

与

（两部分都要计算）；（2）

2 sin?与

cos2?

(1)>>syms x y
>>f1=sqrt(8-x.^2);
>>f2=(1/2)*x.^2;
>>f3=int(f1,-2,2);
>>f4=int(f2,-2,2);
>>f=vpa(f3)-vpa(f4)
f =7.616518051291981025862009923

6．计算半立方抛物线	y	2	?	2 (	x	?	1)	3	被抛物线	y ?	x	截得的一段弧的长度。
				3							3

实验四多元函数微积分

求多元函数的偏导数与全微分

1.1 设

sin(

)

cos

2 xy

求

?x?y

>>clear
>>syms x y
>>f=sin(x*y)+(cos(x*y)).^2;
>>v=[x,y];
>>jacobian(f,v)
ans=[ y*cos(x*y) - 2*y*cos(x*y)*sin(x*y), x*cos(x*y) -2*x*cos(x*y)*sin(x*y)]
>>diff(f,x,2)
ans =- 2*y^2*cos(x*y)^2 + 2*y^2*sin(x*y)^2 - y^2*sin(x*y)
>>diff(diff(f,y),x)
ans=- 2*x*y*cos(x*y)^2 - 2*cos(x*y)*sin(x*y) + cos(x*y) +2*x*y*sin(x*y)^2 - x*y*sin(x*y)

1.2 设

sin

,求

>>syms u v
>>x=exp(u)+u*sin(v);
>>y=exp(u)-u*cos(v);
>>diff(x,u)
ans=exp(u) + sin(v)
>>diff(y,u)
ans=exp(u) - cos(v)
>>diff(x,v)
ans=u*cos(v)
>>diff(y,v)
ans=u*sin(v)
微分学的几何应用

1.3 求出曲面

在点(1,1)处的切平面、法线方程, 并画出图形.

[x,y]=meshgrid(-5:0.1:5);
z=2.*x.^2+y.^2;
mesh(x,y,z)
holdon
[x,y]=meshgrid(-10:0.1:10);
z=4*x+2*y-3;
plot3(x,y,z)
holdon
line([41,-39],[21,-19],[-7,13])
axis([-2020 -20 20 -40 40])

1.4 求曲面	k	(	x	,		)	?				4				在点	??	1	,	1	,		? ?	处的切平面方程, 并把曲面和它的切平面作在同一图形里.
								x	2	?	y	2	?	1		?	4		2		21	?

>>syms x y k;
df_dx=diff(4/(x^2+y^2+1),x)
df_dy=diff(4/(x^2+y^2+1),y)
a=linspace(-10,10,100);

b=a;
[a,b]=meshgrid(a,b);
c=4./(a.^2+b.^2+1);
d=-8/((1/4)^2+(1/2)^2+1)^2*(1/4);
e=-8/((1/4)^2+(1/2)^2+1)^2*(1/2);
f=d.*(a-1/4)+e.*(b-1/2)+/21;
mesh(a,b,c);
holdon;
mesh(a,b,f);
axis([-10,10,-10,10,-2,5]);
df_dx=-(8*x)/(x^2 + y^2 + 1)^2
df_dy=-(8*y)/(x^2 + y^2 + 1)^2

多元函数的极值

1.5 求

(

)

的极值.

>>syms x y;
f=x^3-y^3+3*x^2+3*y^2-9*x;
fx=diff(f,x)
fy=diff(f,y)
fxx=diff(fx,x)
fxy=diff(fx,y)
fyy=diff(fy,y)
fx=3*x^2 + 6*x - 9
fy=6*y - 3*y^2
fxx=6*x + 6
fxy=0
fyy=6 - 6*y
(-3,2)(1,0)

1.6 求函数

在条件

下的极值.

1.6求函数在条件下的极值.

>>syms x y m;
z=x^2+y^2;
df_dy=diff(z,y);
df_dx=diff(z,x);
q=x^2+y^2+x+y-1;
dq_dx=diff(q,x);
dq_dy=diff(q,y);
[x,y,m]=solve(df_dx+m*dq_dx,df_dy+m*dq_dy,q)

3^(1/2)/3- 1
- 3^(1/2)/3 - 1

3^(1/2)/2- 1/2
- 3^(1/2)/2 - 1/2

3^(1/2)/2- 1/2
- 3^(1/2)/2 - 1/2
实验2多元函数积分学（基础实验）
计算重积分

2.1 计算

??D

2dxdy

其中

为由

2 ,

所围成的有界区域.

>>syms x y;int(int(x*y^2,x,2-y,y^0.5),y,1,2)

ans=193/120

2.2 计算

????

(

)

dxdydz

, 其中

由曲面

与

围成.

>>syms t r z;int(int(int((r^2+z)*r,z,r,(2-r^2)^0.5),r,0,1),t,0,2*pi)ans =(pi*(32*2^(1/2) - 25))/30
重积分的应用

2.3 求由曲面

与

所围成的空间区域

的体积.

>>syms t r;
int(int((3/2-r^2)*r,r,0,(3/2)^0.5),t,0,2*pi)
ans=(9*pi)/8

2.4 在

Oxz

平面内有一个半径为2 的圆, 它与

轴在原点

相切, 求它绕

轴旋转一周所得旋转体体积.

>>syms x;
int(4*pi*x*(4-(x-2)^2)^0.5,x,0,4)
ans=16*pi^2
计算曲线积分

2.5 求

(

)

, 其中

积分路径为

3 t

2 .

z t

, 将曲线积分化为定积分)

(注意到,弧长微元

2 t

>>syms t;
x=t;y=t^2;z=3*t^2;f=diff([x,y,z],t);
fun=inline('((1+30*t.^2).^0.5+10*t.^2).*(1+40*t.^2).^0.5','t');quad(fun,0,2)
ans=348.9428

2.6 求

?L. dr

, 其中

( t

)

cos

sin

, 0

(

2 )

>>syms t;
x=cos(t);
y=sin(t);
int(x*y^6*(-2*sin(t))+3*x*(x*y^5+2)*cos(t),t,0,2*pi)ans =6*pi
计算曲面积分

2.7 计算曲面积分

???

(

)

, 其中

为锥面

被柱面

所截得的有限部分.

(注意到,面积微元

dxdy

, 投影曲线

的极坐标方程为

cos

将曲面积分化作二重积分,并采用极坐标计算重积分.)

>>syms t r;

x=r*cos(t);
y=r*sin(t);
z=r;
int(int((x*y+y*z+x*z)*r*2^0.5,r,0,2*cos(t)),t,-pi/2,pi/2)ans =(*2^(1/2))/15

2.8 计算曲面积分

???

dydz

dzdx

dxdy

其中

为球面

的外侧.

symst s r;
symsa real;
int(int(int(3*r^4*sin(s),r,0,a),s,0,pi),t,0,2*pi)ans =(12*pi*a^5)/5
曲线拟合

3.1 为研究某一化学反应过程中温度							x?( C	)	对产品得率			y	(%)	的影响, 测得数据如下:
x	100	110	120	130	140	150	160		170	180	190
y	45	51	54	61	66	70	74		78	85

试求其拟合曲线.

>>clear
>>x=[100,110,120,130,140,150,160,170,180,190];y=[45,51,54,61,66,70,74,78,85,];
a=polyfit(x,y,1)
z=polyval(a,x);
plot(x,y,'gp',x,z,'r');
a= 0.4830 -2.7394

3.2给定平面上点的坐标如下表:

x	0 . 1	0 . 2	0 . 3	0 . 4	0 . 5	0 . 6	0 . 7	0 . 8	0 . 9
y	5 . 1234	5 . 3057	5 . 5687	5 . 9378	6 . 4337	7 . 0978	7 . 9493	9 . 0253	10 . 3627

试求其拟合曲线.

>>x=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9];
y=[5.1234,5.3057,5.5687,5.9378,6.4337,7.0978,7.9493,9.0253,10.3627];a=polyfit(x,y,3)
z=polyval(a,x);
plot(x,y,'bp',x,z,'r');

a=4.9875 0.6902 1.3202 4.9774

3.3已知lg(x)的[1，101]区间11个整数采样点的函数值如下表

x	1	11	21	31	41	51
lg(x)	0	1.0414	1.3222	1.4914	1.6128	1.7076

试求lg(x)的5次拟合多相式p(x)，并分别绘制出lg(x)和p(x)在[1，101]区间的函数曲线

x=[111 21 31 41 51];

lgx=[01.0414 1.3222 1.4914 1.6128 1.7076];

P=polyfit(x,lgx,5)

X=[1:0.2:51];

Y=log10(X);

plot(X,Y,'*',X,polyval(P,X),'-')

gridminor

实验六无穷级数及微分方程（基础实验）

数项级数

1.1(1)观察级数	? ?	1		的部分和序列的变化趋势.
	n?1	n	2

x=0;
forn=1:50;
x=x+1/n^2;
plot(n,x,'r*');
hold on
end

(2) 观察级数	? ?	1	的部分和序列的变化趋势.
	n?1	n

x=0;
forn=1:100;
x=x+1/n;
plot(n,x,'r*')
hold on
end

1.2 设	a	?	10	n	,	求	? ?	a		.
		n	n !				n?1		n

s=10;

fori=1:inf;

s=s+s*10/(i+1);

end

s=5.2257e+086

求幂级数的收敛域

1.3 求

?
?

(

3 )

的收敛域与和函数.

>>syms n x k;

limit(4^(2*n)/(n+1)/(4^(2*n+2)/(n+2)),n,inf)

s=symsum(4^(2*n)*(x-3)^n/(n+1),n,0,inf)

ans=1/16

s=piecewise([49/16 <= x, Inf], [x <> 49/16 and abs(16*x - 48)<= 1, -log(49 - 16*x)/(16*x - 48)])

1.4 设	a	?	3	n	,	求	? ?	a		.(*)
		n	n !				n?1		n

symsn;S=symsum(3^sym('n!'),n,1,inf)

1.5求下列级数的和：

?
?

，

?
?

，

?
?

(

?1 )

，

100
?

n?1

2 (

1 )

n?1

(1)

>>syms n

>>symsum((2*n-1)/2.^n,1,inf)

ans=3

(2)

>>syms n

>>symsum(1/(2*(2*n+1)),1,inf)

ans=Inf
(3)
>>syms n
>>symsum((-1)^(n+1)/n,1,inf)
ans=log(2)
(4)
>>syms n
>>symsum(n.^2,1,100)
ans =338350

求幂级数的收敛域

1.6 求	? ?	(	?1 )	n	?1	x	n	的收敛域与和函数.
	n?1					n

>>syms n x.

>>s=symsum((-1)^(n-1)*(x.^n/n),n,1,inf)
s=piecewise([x = -1, -Inf], [abs(x) <= 1 and x <> -1, log(x +1)])函数的幂级数展开

1.7将函数sinx展开为x的幂级数，分别展开至5次和20次。

>>syms x
>>f=sin(x);
>>f1=taylor(f,x,5)
f1=x - x^3/6
>>f2=taylor(f,x,20)
f2=- x^19/1215100408832000 + x^17/355687428096000 -x^15/1307674368000 + x^13/6227020800 - x^11/39916800 + x^9/362880 -x^7/5040 + x^5/120 - x^3/6 + x

1.8 将函数(1

展开为x的幂级数，m为任意常数。展开至4 次幂。

>>syms m x
>>f=(1+x).^m;
>>f1=taylor(f,x,4)
f1=(m/3 - m*(m/4 - m^2/6) - m^2/4)*x^3 + (m^2/2 - m/2)*x^2 + m*x + 1

f x ( )

1.9 将函数

展开为(

x ?

的幂级数。

>>syms x

>>f=1/(x.^2+5*x-3);

>>f1=taylor(f,x-2)

f1=taylor(1/(x^2 + 5*x - 3), x = x - 2, 6, AbsoluteOrder)

1.10 将函数cos x展开成

(

?
)

的幂级数，取前10 项。

symsx;

f=cos(x);
taylor(f,x,10,-pi/3)
ans=
1/2+1/2*3^(1/2)*(x+1/3*pi)-1/4*(x+1/3*pi)^2-1/12*3^(1/2)*(x+1/3*pi)^3+1/48*(x+1/3*pi)^4+1/240*3^(1/2)*(x+1/3*pi)^5-1/1440*(x+1/3*pi)^6-1/10080*3^(1/2)*(x+1/3*pi)^7+1/800*(x+1/3*pi)^8+1/725760*3^(1/2)*(x+1/3*pi)^9

1.11 求函数

f x ( )

在[

???

上的傅立叶级数。

1.12 求出函数

f x ( )

在区间[

???

上的前11 个傅立叶系数，即n＝5。

function
[a0,ak,bk]=myfly(f);syms k x
a0=int(f,x,-pi,pi)/pi;ak=int(f*cos(k*x),x,-pi,pi)/pi;
bk=int(f*sin(k*x),x,-pi,pi)/pi;
%fourieran.m
functionan=fourieran(f,n) syms x an=int(f*cos(n*x),x,-pi,pi)/pi;
%fourierbn.m
functionbn=fourierbn(f,n) syms x bn=int(f*sin(n*x),x,-pi,pi)/pi; clear syms xn f=x^2+x^3 a0=fourieran(f,0); for n=1:5
a(n)=fourieran(f,n); end for n=1:5
b(n)=fourierbn(f,n); end
f= x^3 + x^2
>> a
a=
[-4, 1, -4/9, 1/4, -4/25]
f =
x^3+ x^2
>>b b =

[2*pi^2 - 12, ((3*pi)/2 - pi^3)/pi, (2*pi^2)/3 - 4/9, 3/16 - pi^2/2,(2*pi^2)/5 - 12/125]

1.13 求

arctan

的5 阶泰勒展开式.

symsx
f=atan(x);
f1=taylor(f,x,5)

f1=

x/26+ atan(5) - (5*(x - 5)^2)/676 + (37*(x - 5)^3)/263 - (15*(x -5)^4)/57122 + (719*(x - 5)^5)/14851720 - 5/26

1.14 求

?1??

?1?2

在

处的8 阶泰勒展开, 并通过作图比较函数和它的近似多

>>clear
>>syms x
>>f=exp(-(x-1).^2*(x+1).^3);
>>f1=taylor(f,1,8)
f1=1/exp(35721)

1.15 利用泰勒公式

R n

(

)

近似计算

. 若

1 ,

要求截断误差

0 . 005 ,

问n 应取多大?

2 !

n !

1.16 观察函数

(

)

sin

各阶泰勒展开的图形.

(1) 固定

x 0?

0 ,

观察阶数n 的影响;

(2) 扩大显示区间范围, 以观察在偏离展开点

x 0

时泰勒多项式对函数的逼近情况;

(3) 固定

10 ,

观察

x 0

的影响.

（1）>>syms x
>>f=sin(x);
>>f1=taylor(f,0)
f1 =x^5/120 - x^3/6 + x
>>ezplot(f1)
>>hold on
>>f2=taylor(f,10)
f2=x^9/362880 - x^7/5040 + x^5/120 - x^3/6 + x
>>ezplot(f2)

求解微分方程

2.1 求微分方程

的通解.

>>y=dsolve('Dy+2*x*y=x*exp(-x^2)','x')

y=C2/exp(x^2) + x^2/(2*exp(x^2))

2.2 求微分方程

x y

在初始条件

y x

?1?

2 e

下的特解.

>>y=dsolve('x*Dy+y=exp(x)','y(1)=2*exp(1)','x') y =(exp(1) + exp(x))/x

2.3 求解微分方程

, 并作出其积分曲线.

>>y=dsolve('D2y-2*x=exp(x)','x')
y =C7 + exp(x) + C6*x + x^3/3
积分曲线
>>x=-2:0.1:2;
>>y=1./3*x.^3+exp(x)+x+1;
>>plot(x,y)

2.4 求微分方程组

???????

dx
?x

dy
?

在初始条件

x t

1 ,

y t

下的特解.

>>[x y]=dsolve('Dx+x+2*y-exp(t)','Dy-x-y','x(0)=1','y(0)=0','t') x=cos(t)
y=exp(t)/2 - cos(t)/2 + sin(t)/2
2.5求出初值问题

sin

cos

(

0 )

1 ,

0 )

的数值解,并作出数值解的图形.

functiondy=ffer(x,y);
dy=zeros(2,1);dy(1)=y(2);
dy(2)=-y(2)*(sin(x))^2-y(1)+(cos(x))^2;
[X,Y]=ode23s('ffer',[04],[1 0]);
plot(X,Y(:,1),'-');

X=
0
0.0183
0.0366

0.0550
0.0733
0.0917
0.1100
0.1284
0.1469
0.1657
0.1870
0.2111
0.2384
0.2694
0.3048
0.3453
0.3920
0.4461
0.5098
0.5863
0.6823
0.8161
0.9767
1.1374
1.2859
1.4337
1.5851
1.7415
1.9018
2.0632
2.2229
2.3794
2.5324
2.6836
2.8377
3.0051
3.1122
3.2193
3.3242
3.4361
3.5540
3.6787
3.8127
3.9528
4.0000
Y=
1.0000 0

1.0000 -0.0000
1.0000 -0.0000
1.0000 -0.0001
1.0000 -0.0001
1.0000 -0.0003
1.0000 -0.0004
1.0000 -0.0007
1.0000 -0.0010
0.9999 -0.0015
0.9999 -0.0022
0.9998 -0.0031
0.9997 -0.0044
0.9996 -0.00
0.9993 -0.0092
0.9988 -0.0132
0.9981 -0.0191
0.9968 -0.0277
0.9947 -0.0405
0.9908 -0.0597
0.9837 -0.06
0.9683 -0.1414
0.9398 -0.2147
0.90 -0.2928
0.8503 -0.3618
0.7924 -0.4213
0.7249 -0.4684
0.92 -0.4994
0.5680 -0.5110
0.4862 -0.5015
0.4084 -0.4711
0.3384 -0.4210
0.2791 -0.3527
0.2321 -0.2671
0.19 -0.1628
0.1823 -0.0343
0.1833 0.0525
0.1936 0.1394
0.2125 0.2209
0.2417 0.2995
0.2811 0.3678
0.3302 0.4181
0.3881 0.4421
0.4495 0.4316
0.4696 0.4202

2.6洛伦兹(Lorenz)方程组是由三个一阶微分方程组成的方程组.这三个方程看似简单,也没有包含复杂的函数,但它的解却很有趣和耐人寻味.试求解洛伦兹方程组

???????

?( t

)

( t

)

( t

)

( t

)
,

?( t

)

( t

)

( t

)

( t

?( t

)

( t

)

( t

)

( t

)

(

0 )

12 ,

(

0 )

4 ,

(

0 )

并画出解曲线的图形.

functiondy=lorenz(t,y)
dy=zeros(3,1);dy(1)=16*y(2)-16*y(1);
dy(2)=-y(1)*y(3)+45*y(1)-y(2);dy(3)=y(1)*y(2)-4*y(3);

[T,Y]=ode45('lorenz',[00.1],[12 4 0])
plot(T,Y(:,1),'-',T,Y(:,2),'*',T,Y(:,3),'+')
T=
0
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0000
0.0001
0.0001
0.0001
0.0001
0.0003
0.0004
0.0005
0.0007

0.0013
0.0020
0.0026
0.0033
0.0058
0.0083
0.0108
0.0133
0.0158
0.0183
0.0208
0.0233
0.0258
0.0283
0.0308
0.0333
0.0358
0.0383
0.0408
0.0433
0.0458
0.0483
0.0508
0.0533
0.0558
0.0583
0.0608
0.0633
0.0658
0.0683
0.0708
0.0733
0.0758
0.0783
0.0808
0.0833
0.0858
0.0883
0.0908
0.0933
0.0950
0.0966
0.0983
0.1000

Y=
12.0000 4.0000 0
11.9999 4.0006 0.0001
11.9997 4.0011 0.0001
11.9996 4.0017 0.0002
11.9995 4.0022 0.0002
11.9988 4.0050 0.0005
11.9981 4.0079 0.0007
11.9975 4.0107 0.0010
11.9968 4.0135 0.0012
11.9934 4.0275 0.0025
11.9901 4.0415 0.0037
11.9868 4.0555 0.0050
11.9835 4.0695 0.0063
11.9670 4.1395 0.0127
11.9507 4.2093 0.0192
11.9346 4.2790 0.0258
11.9187 4.3486 0.0326
11.8416 4.6950 0.0676
11.7690 5.0386 0.1049
11.7007 5.3797 0.1446
11.6367 5.7182 0.1865
11.4298 6.9914 0.3669
11.2804 8.2360 0.5783
11.1853 9.4575 0.8197
11.1413 10.6612 1.0908
11.1460 11.8517 1.3920
11.1970 13.0336 1.7238
11.2923 14.2107 2.0877
11.4299 15.3865 2.4850
11.6083 16.54 2.9180
11.8260 17.7471 3.30
12.0816 18.9369 3.9006
12.3740 20.1357 4.4560
12.7022 21.3451 5.0584
13.0651 22.5660 5.7117
13.4620 23.7988 6.4197
13.18 25.0437 7.1867
14.3538 26.2998 8.0170
14.8472 27.5659 8.9154
15.3710 28.8399 9.8867
15.9244 30.1193 10.9358
16.5063 31.4003 12.0678
17.1155 32.6786 13.2876

17.7509 33.9488 14.6002
18.4108 35.2047 16.0103
19.0938 36.4388 17.5221
19.7977 37.28 19.1398
20.5206 38.8071 20.8666
21.2598 39.9212 22.7049
22.0125 40.9734 24.6563
22.7755 41.9510 26.7213
23.5453 42.8405 28.86
24.3177 43.6275 31.1856
25.0885 44.2970 33.5776
25.8527 44.8337 36.0682
26.6052 45.2222 38.83
27.3402 45.4476 41.3068
27.8220 45.5001 43.1334
28.2913 45.4679 44.9845
28.7462 45.3473 46.8548
29.1845 45.1348 48.7385

实验七矩阵运算与方程组求解

1 计算范德蒙行列式	1	1		1	1		1	.
	x 1	x	2	x 3	x	4	x 5
	x 1 2	x	2 2	x 3 2	x	2 4	x 5 2
	x 1 3	x	3 2	x 3 3	x	3 4	x 5 3
	x 1 4	x	4 2	x 3 4	x	4 4	x 5 4

>>syms x1 x2 x3 x4 x5
>>
A=[1,1,1,1,1;x1,x2,x3,x4,x5;x1^2,x2^2,x3^2,x4^2,x5^2;x1^3,x2^3,x3^3,x4^3,x5^3;x1^4,x2^4,x3^4,x4^4,x5^4];
>>a=det(A);
>>a=simple(a)
a =(x1 - x2)*(x1 - x3)*(x1 - x4)*(x2 - x3)*(x1 - x5)*(x2 - x4)*(x2 -x5)*(x3 - x4)*(x3 - x5)*(x4 - x5)

2 设矩阵	A	?	????????	3	7	2		6	?4		? ? ? ??, ? ? ?	求	\|	A	\|,	tr	(	A ),	A	3	.
				7	9	4		2	0
				11	5	?	6	9	3
				2	7	?	8	3	7
				5	7	9		0	?	6

>>A=[3 7 2 6 -4;7 9 4 2 0;11 5 -6 9 3;2 7 -8 3 7;5 7 9 0 -6]; >>det(A)
ans=11592
>>A'
ans=

3 7 11 2 5
7 9 5 7 7
2 4 -6 -8 9
6 2 9 3 0
-4 0 3 7 -6
>>A^3
ans=
Columns 1 through 4

726 2062 944 294
1848 3150 26 1516
1713 2218 31 1006
1743 984 -451 1222
801 2666 477 745
Column5

-358
228
404
384
-125

3 设

?????

求矩阵M 的秩.

>>m=[3 2 -1 -3 -2;2 -1 3 1 -3;7 0 5 -1 -8];
rank(m)
ans= 2

4 已知矩阵	M	?	?????	3	2		?	1	?	3	?????	的秩等于2, 求常数t 的值.
				2	?	1	3		1
				7	0		t		?	1

symst;
M=[32 -1 -3 1;2 -1 3 1 1;7 0 t -1 1]
m=rref(M) %分母为t-5，将5代入M

M=[32 -1 -3 1;2 -1 3 1 1;7 0 5 -1 1];
refm=rref(M)

[3, 2, -1, -3, 1]
[2, -1, 3, 1, 1]
[7, 0, t, -1, 1]

[1, 0, 0, -1/7, (3*t - 5)/(7*t - 35)]
[0, 1, 0, -9/7, -(t + 17)/(7*t - 35)]
[0, 0, 1, 0, -2/(t - 5)]

refm=

Columns1 through 4

1.0000 0 0.7143 -0.1429
0 1.0000 -1.5714 -1.2857
0 0 0 0

Column5

0
0
1.0000
t=5

5 设

?????

?
?

求矩阵A 的秩.

>>A=[2 -3 8 2;2 12 -2 12;1 3 1 4];
rank(A)
ans=2

6 求向量组

( 1 ,

2 ,

?1 , 1 ),?3

(

0 ,

?4 , 5 ,

?2 ),?2

(

2 , 0 , 3 , 0 )

的秩.

>>a1=[1 2 -1 1];
a2=[0-4 5 -2];
a3=[20 3 0];

rank([a1;a2;a3])
ans= 2
7求向量组

( 1 ,

?1 ,

2 ,

4 ),?2

(

0 , 3 , 1 ,

2 ),?3

( 3 , 0 , 7 , 14 ),?4

( 1 ,

?1 ,

2 , 0 ),?5

(

2 , 1 , 5 , 0 )

的极大无关组,并将其它向量用极大无关组线性表示.

>>a1=[1 -1 2 4];
a2=[03 1 2];
a3=[30 7 14];
a4=[1-1 2 0];
a5=[21 5 0];
rank([a1;a2;a3;a4;a5])
rank([a1;a2;a3])
rank([a1;a3;a4])
rank([a1;a2;a4])
ans=3
ans=2
ans=3
ans=3

8 求解线性方程组	???????	x 1			?		x	2	?		2	x	3	?		x	4		?		0 ,
		3		x	1	?	x	2	?		x	3	?		2	x	4		?		0 ,
		5			x	2	?		7	x	3	?		3	x	4	?			0 ,
		2	x 1		?		3	x	2	?		5	x	3	?		x	4		?		0 .

>>a=[1 1 -2 -1;3 -1 -1 2;0 5 7 3;2 -3 -5 -1]; b=[0 0 0 0]';
>>x=a\b

0
0
0
0

9 向量组

( 1 , 1 ,

2 , 3 ),?2

( 1 ,

?1 , 1 , 1 ),?3

( 1 , 3 ,

4 , 5 ),?4

( 3 , 1 , 5 , 7

)

是否线性相关?

>>a1=[1 1 2 3];
a2=[1-1 1 1];
a3=[13 4 5];
a4=[31 5 7];
rank([a1;a2;a3;a4])
ans=3
线性无关

10 求出通过平面上三点(0,7),(1,6)和(2,9)的二次多项式

并画出其图形.

>>A=[0 0 1;1 1 1;4 2 1];

>>B=[7 6 9]';
>>abc=inv(A)*B

abc=
2
-3
7
>>ezplot('2*x^2-3*x+7')

11 求出通过平面上三点(0,0),(1,1),(-1,3)以及满足

?1 )

20 ,

?( 1 )

的4 次多项式

12 解方程组

???????

x 1

>>a=[1 -1 2 1;2 -1 1 2;1 0 -1 1;3 -1 0 3];
b=[1;3;2;5];
rref([ab])

ans=

1 0 -1 1 2
0 1 -3 0 1
0 0 0 0 0
0 0 0 0 0
x1=2+x3-x4;

x2=1+3*x3;

13 当

为何值时,方程组

?????

ax 1

无解、有唯一解、有无穷多解?当方程组有

x 1

解时,求通解.

>>syms a
>>A=[a 1 1;1 a 1;1 1 a];

>>Ab=[a 1 1 1;1 a 1 1;1 1 a 1];

>>b=[1 1 1]';

>>rref(A)

ans=

[1, 0, 0]

[0, 1, 0]

[0, 0, 1]

>>rref(Ab)

ans=

[1, 0, 0, 1/(a + 2)]

[0, 1, 0, 1/(a + 2)]

[0, 0, 1, 1/(a + 2)]

>>Ab1=[-2 1 1 1;1 -2 1 1;1 1 -2 1];

>>rref(Ab1)

ans=

1 0 -1 0

0 1 -1 0

0 0 0 1

>>x=inv(A)*b

(a+ 1)/(a^2 + a - 2) - 2/(a^2 + a - 2)

实验八矩阵的特征值与特征向量

1 求矩阵

?????

?1?.

的特征值与特值向量.

>>A=[-1 0 2;1 2 -1;1 3 0];

>>[V D]=eig(A,'nobalance')

Columns1 through 2

1.0000 1.0000 - 0.0000i
-0.3333 -0.0000 + 0.0000i
-0.0000 1.0000 - 0.0000i

Column3

1.0000+ 0.0000i
-0.0000 - 0.0000i
1.0000 + 0.0000i

Columns1 through 2

-1.0000 0
0 1.0000 + 0.0000i
0 0

Column3

0
0
1.0000 - 0.0000i

2 已知

( 1 , 1 ,

?1 )

是方阵

?????

的一个特征向量,求参数

及特征向量

所属的特征值.

>>clear
>>syms a b c;
>>A=[2-c -1 2;5 a-c 3;-1 b -2-c];
>>x=[1;1;-1];
>>A*x
ans=
- c - 1
a - c + 2
b + c + 1
>>[a b c]=solve('-1-c','2+a-c','1+b+c','a,b,c') a =-3
b=0
c=-1

3 设矩阵	A	?	?????	4	1	1	?????	,求一可逆矩阵	P	,使	P1	AP	为对角矩阵.
				2	2	2
				2	2	2

>>A=[4 1 1;2 2 2;2 2 2];
>>[P D]=eig(A)

0.5774 0.5774 0
0.5774 -0.5774 -0.7071
0.5774 -0.5774 0.7071

6.0000 0 0
0 2.0000 0
0 0 0

4 已知方阵

??2

?????

与

??1

?????

相似, 求

>>syms x;
>>A=[-2 0 0;2 x 2;3 1 1];
>>E=eig(A)
E=
-2
x/2 - (x^2 - 2*x + 9)^(1/2)/2 + 1/2
x/2 + (x^2 - 2*x + 9)^(1/2)/2 + 1/2
>>y=-2
y= -2
>>x=solve(1/2*x+1/2+1/2*(x^2-2*x+9)^(1/2)-2) x =0

5 求一个正交变换,化二次型

x 1

x 3

为标准型.

>>A=[0 1 1 0;1 0 1 0;1 1 0 0;0 0 0 2];
>>[Q D]=eig(A)

-0.7152 0.3938 0.5774 0
0.0166 -0.8163 0.5774 0
0.6987 0.4225 0.5774 0
0 0 0 1.0000

-1.0000 0 0 0
0 -1.0000 0 0
0 0 2.0000 0
0 0 0 2.0000

6已知二次型

f	(	x 1	,	x	2	,	x	3	)	?	x 1 2	?	2	x	2	?	x 3 2	?	2	x 1	x	2	?	4	x 1	x	3	?	2	x	2	x	3
															2

(1)求标准形; (2)求正惯性指数; (3)判断二次型是否正定.

（1）
>>A=[1 1 -2;1 -2 1;-2 1 1];
>>[Q D]=eig(A)

0.4082 0.5774 -0.7071
-0.8165 0.5774 -0.0000
0.4082 0.5774 0.7071

-3.0000 0 0
0 0.0000 0
0 0 3.0000
（2）由对角矩阵D得，正惯性指数是1。

（3）>>D=diag([-3,0,3]);
ifall(D>0)
disp('二次型正定')
elsedisp('二次型非正定')
end
二次型非正定

中南大学 数学建模报告

中南大学数学建模报告