数学模型习题解答解读

上机练习题一

班级： 姓名： 学号：

1.建立起始值=3,增量值=5.5,终止值=44的一维数组x 答案： x=(3:5.5:44)

2.写出计算 Sin(30o )的程序语句. 答案： sin(pi*30/180) 或 sin(pi/6)

3.矩阵??????????=187624323A ,矩阵????

??????=333222111B ;分别求出B A ?及A 与B 中对应元素之间的乘积. 答案：A = [3,2,3; 4,2,6; 7,8,1]

B = [1,1,1; 2,2,2; 3,3,3]

A*B ；A.*B

4计算行列式的值1

876243

23=A 。答案：det(A)

5对矩阵 ????

??????=187624323A 进行下述操作。 (1)求秩。答案：rank(A)

(2)求转置。答案：A'

(3) 对矩阵求逆，求伪逆。答案：inv(A) ，pinv(A)

(4) 左右反转,上下反转。答案：fliplr(A)，flipud(A)

(5) 求矩阵的特征值. 答案：[u,v]=eig(A)

(6) 取出上三角和下三角. 答案：triu(A) tril(A)

(7)以A 为分块作一个3行2列的分块矩阵。答案：repmat(a)

6 计算矩阵??????????897473535与????

??????638976242之和。

>> a=[5 3 5;3 7 4;7 9 8];

>> b=[2 4 2;6 7 9;8 3 6];

>> a+b

7 计算??????=572396a 与??

????=864142b 的数组乘积。 >> a=[6 9 3;2 7 5]; >> b=[2 4 1;4 6 8];

ans =

12 36 3

8 42 40

8 已知：????

??????=987654321a ，分别计算a 的数组平方和矩阵平方，并观察其结果。

>> a=[1 2 3;4 5 6;7 8 9];

>> a.^2

ans =

1 4 9

16 25 36

49 64 81

>> a^2

ans =

30 36 42

66 81 96

102 126 150

上机练习题二

班级： 姓名：

学号：

1 对于B AX =，如果??????????

=753467294A ，????

?

?????=282637B ，求解X 。

>> A=[4 9 2;7 6 4;3 5 7];

>> B=[37 26 28]’;

>> X=A\B

X =

-0.5118

4.0427

1.3318

2 角度[]604530=x ，求x 的正弦、余弦、正切和余切。

>> x=[30 45 60];

>> x1=x/180*pi;

>> sin(x1)

ans =

0.5000 0.7071 0.8660

>> cos(x1)

ans =

0.8660 0.7071 0.5000 >> tan(x1)

0.5774 1.0000 1.7321

>> cot(x1)

ans =

1.7321 1.0000 0.5774

3 将矩阵??

????=7524a 、??????=3817b 和??????=2695c 组合成两个新矩阵： （1）组合成一个4?3的矩阵，第一列为按列顺序排列的a 矩阵元素，第二列为按列顺序排列的b 矩阵元素，第三列为按列顺序排列的c 矩阵元素，即

????

?

?

??????237912685574

（2）按照a 、b 、c 的列顺序组合成一个行矢量，即

[]296531877254

答案：a=[4 2;5 7];

b=[7 1;8 3];

c=[5 9;6 2];

% (1)

>> d=[a(:) b(:) c(:)]

d =

4 7 5

5 8 6

2 1 9

7 3 2

% (2)

>> e=[a(:);b(:);c(:)]'

e =

4 5 2 7 7 8 1 3 5 6

9 2 或利用（1）中产生的d

>> e=reshape(d,1,12)

ans =

4 5 2 7 7 8 1 3 5 6

9 2 4求解在x =8时多项式(x -1)(x -2) (x -3)(x -4)的值。

>> p=poly([1 2 3 4]);

>> polyvalm(p,8)% polyval(p,8)

ans =

840

5求方程023973234=-++x x x 的全部根。

p=[3,7,9,0,-23]; %建立多项式系数向量

x=roots(p) %求根

上机练习题三

班级： 姓名： 学号：

1、 设x 是数组，求均值和方差

解：函数文件如下：

function [xx,s]=func1(x)

n=length(x);

xx=sum(x)/n;

s=sqrt((sum(x.^2)-n*xx^2)/(n-1));

命令窗口：

>> x=[1 2 3 4 5];[xx,s]=func1(x)

2、求满足100)1ln(0

>+∑=m

n n 的最小m 值

s=0;

n=0;

while(s<=100)

s=s+log(1+n);

n=n+1;

end

n,s

3、用循环语句形成Fibonacci 数列,....4,3,,12121=+===--k F F F F F k k k 。并验证极限2

5

11+→

-k k F F （

提示：计算至两边误差小于精度1e-8为止）

解： 求Fibonacci 数列的函数文件：

function fibo(k)

clc;

f(1)=1;

f(2)=1;

for n=3:k

f(n)=f(n-1)+f(n-2)

end

f(k)

验证极限的函数文件：

function [k,a]=funTest(e)

a=abs(1-(1+sqrt(5))/2);

k=2; while(a>e)

k=k+1; a=abs(fun(k)/fun(k-1)-(1+sqrt(5))/2); end 命令行：

>> [k,a]=funTest(10^-8)

k =

21

a =

9.7719e-009

或者M 文件如下：

clear; F(1)=1;F(2)=1;k=2;x=0;

e=1e-8; a=(1+sqrt(5))/2;

while abs(x-a)>e

k=k+1; F(k)=F(k-1)+F(k-2); x=F(k)/F(k-1);

end

a,x,k 4、分别用for 和while 循环结构编写程序，求出∑==

610123i i

K ，并考虑一种避免循环语句的程序设计，比较各种算法的运行时间。

解：for 循环结构：M 文件loop.m

k=0;

for i=1:10^6

k=k+sqrt(3)*2^-i;

end

k

while 循环结构：M 文件loop1.m

k=0;i=1;

while i<=10^6

k=k+sqrt(3)*2^(-i);

i=i+1;

end

k

非循环结构：M 文件nonLoop.m

i=1:10^6;

x=sqrt(3)*(2.^-i);

k=sum(x) 速度比较：>>tic;loop;toc ％循环结构的执行时间

k =

1.7321

Elapsed time is 1.813000 seconds.

>> tic;nonLoop;toc ％非循环结构的执行时间

k =

1.7321

Elapsed time is 1.094000 seconds.

上机练习题四

班级： 姓名： 学号：

1、作图描述气温变化

>> x=0:24;

>> y=[15,14,14,14,14,15,16,18,20,22,23,25,28,31,32,31,29,27,25,24,22,20,18,17,16];

>> plot(x,y)

2、作出下列函数图形

（1）)2sin(22--=x x x y 22≤≤-x （分别使用plot 和fplot 完成）

解：>> fplot('x^2*sin(x^2-x-2)',[-2 2]) ％fplot 方法

>> x=-2:0.1:2;y=x.^2.*sin(x.^2-x-2);plot(x,y) ％plot 方法 如图（4.1）

（2）19

42

2=+y x （椭圆 提示：用参数方程） 解：>> r=-pi:0.1:pi;x=2*cos(r);y=3*sin(r);plot(x,y) ％ 如图（4.2）

解法二

x=-2:1/100:2;

y1=3*sqrt(1-x.^2/4); y2=-3*sqrt(1-x.^2/4);

plot(x,y1,'r-',x,y2,'r-'); axis equal tight;

图（4.1） 图（4.2）

（3）22y x z += （抛物面） 3,3<

解：（错误）>> x=[-3:0.1:3];y=[-3:0.1:3];z=x.^2+y.^2; plot3(x,y,z) ％ 如图（4.31）

（正确）>> xa=-3:0.1:3;ya=-3:0.1:3;[x,y]=meshgrid(xa,ya); ％ 如图（4.32）

>> z=x.^2+y.^2;mesh(x,y,z); >> surf(x,y,z)

图（4.31）error 图（4.32）

（4）曲面133,3,622232224<<-<+---++=y x y x y x y x x z

解： >> xa=linspace(-3,3,100);ya=linspace(-3,13,100);

>> [x,y]=meshgrid(xa,ya);

>> z=x.^4+3*x.^2+y.^2-2*x-2*y-2*x.^2.*y+6;

>> mesh(x,y,z)

>> surf(x,y,z)

（5）空间曲线20),2cos(,cos ,sin <<===t t z t y t x

解：>> t=linspace(0,2,50);x=sin(t);y=cos(t);z=cos(2*t);

>> plot3(x,y,z)

（6）半球面2sin cos ,2sin sin ,2cos ,0360,090x y z φθφθφθφ===≤≤≤≤

解： >> a=linspace(0,pi/2,50);b=linspace(0,2*pi,50);

>> [a,b]=meshgrid(a,b); >> x=2*sin(a).*cos(b);y=2*sin(a).*sin(b);z=2*cos(a);

>> surf(x,y,z)

（7）三条曲线合成图π<<-===x x y x x y x y 0,sin ),10sin(sin ,sin 321

解： >>x=linspace(0,pi,50);y1=sin(x);

>> plot(x,y1);hold on;

>> y2=sin(x).*sin(10*x);

>> plot(x,y2);

>> y3=-sin(x);

>> plot(x,y3);

>> hold off;

3、作下列分段函数图??

???-<≤>-=1.11.11.11.11.1x x x x y

x=-5:0.1:5;

for i=1:length(x)

if x(i)>1.1

y(i)=1.1;

elseif x(i)<-1.1

y(i)=-1.1;

else

y(i)=x(i);

end

end

plot(x,y);

grid on;

或者M 文件如下：

x=-5:0.1:5;

y=1.1.*( x>1.1)+x.*( x <=1.1&x>=-1.1)-1.1.*(x<-1.1);

plot(x,y)

4、用MATLAB 函数表示下列函数，并作图。

??

???-≤++--≤+<--->+---=1),5.175.375.0exp(5457.011),6exp(7575.01),5.175.375.0exp(5457.0),(222222y x x x y y x x y y x x x y y x p

解：建立M 文件pxy 如下：

xa=-2:0.05:2;ya=xa;

nx=length(xa);ny=length(ya);

[x,y]=meshgrid(xa,ya);

z=zeros(nx,ny);

[a1,b1]=find(x+y>1); %第a1列b1行对应的x+y>1 (x 对应列；y 对应行)

%第a1列对应的x 值是xa(a1)；第b1行对应的y 值是ya(b1)

z((a1-1)*ny+b1)=0.5457*exp(-0.75*ya(b1).^2-3.75*xa(a1).^2-1.5*xa(a1));

[a2,b2]=find(x+y<=1&x+y>-1);

z((a2-1)*ny+b2)=0.7575*exp(-ya(b2).^2-6*xa(a2).^2);

[a3,b3]=find(x+y<=-1);

z((a3-1)*ny+b3)=0.5457*exp(-0.75*ya(b3).^2-3.75*xa(a3).^2+1.5*xa(a3));

surf(x,y,z);

命令窗口：

>> pxy

运行结果如右图：

或者M文件如下：

clear;close;

xa=-2:0.1:2;ya=-2:0.1:2;[x,y]=meshgrid(xa,ya);

z=zeros(size(x));

k1=find(x+y>1);

z(k1)=0.5457*exp(-0.75*y(k1).^2-3.75*x(k1).^2-1.5*x(k1));

k2=find(x+y<=1&x+y>-1);

z(k2)=0.7575*exp(-y(k2).^2-6*x(k2).^2);

k3=find(x+y<-1);

z(k3)=0.5457*exp(-0.75*y(k3).^2-3.75*x(k3).^2+1.5*x(k3)); mesh(x,y,z);

或者M文件如下：

xa=-2:0.05:2;ya=xa;

[x,y]=meshgrid(xa,ya);

z=0.5457*exp(-0.75*y.^2-3.75*x.^2-1.5*x).*(x+y>1)+...

0.7575*exp(-y.^2-6*x.^2).*(x+y<=1&x+y>-1)+...

0.5457*exp(-0.75*y.^2-3.75*x.^2+1.5*x).*(x+y<=-1);

surf(x,y,z);

上机练习题五

班级：姓名：学号：

1、运行demo

解：>>demo

2、查询trapz的功能、用法、目录、程序结构、相同目录下其它文件

解：>> help trapz ――功能用法

>> type trapz――程序结构，源码

>> which trapz――所在目录

>> help C:\MATLAB6p5\toolbox\matlab\datafun――该目录下其它文件

3在[0,4pi]画sin(x),cos(x)(在同一个图象中); 其中cos(x)图象用红色小圆圈画.并在函数图上标注“y=sin(x)”, “y=cos(x)” ,x轴,y轴,标题为“正弦余弦函数图象”.

x=linspace(0,4*pi,100);

y=sin(x);

plot(x,y);

gtext('y = sin(x)');% 图形注解,注意要用鼠标定位

hold on;

y=cos(x);

plot(x,y,'ro');

gtext ('y = cos(x)');% 图形注解

xlabel('x轴'); % x轴注解

ylabel('y轴'); % y轴注解

title('正弦余弦函数图象'); % 图形标题

4从键盘输入若干个数，当输入0时结束输入，求这些数的平均值和它们之和。

sum=0;

cnt=0;

val=input('Enter a number (end in 0):');

while (val~=0)

sum=sum+val;

cnt=cnt+1;

val=input('Enter a number (end in 0):');

end

if (cnt > 0)

sum

mean=sum/cnt

end

5若一个数等于它的各个真因子之和，则称该数为完数，如6=1+2+3，所以6是完数。求[1,500]之间的全部完数。

for m=1:500

s=0;

for k=1:m/2

if rem(m,k)==0

s=s+k;

end

end

if m==s

disp(m);

end

end

上机练习题六

班级： 姓名： 学号：

1 假定数据点来源为：x e x x x f x sin )53()(52-+-=，试根据生成的数据进行插值处理，得出较平滑的曲线。 x=0:.12:1;y=(x.^2-3*x+5).*exp(-5*x).*sin(x);plot(x,y,x,y,'o')

x1=0:.02:1;y0=(x1.^2-3*x1+5).*exp(-5*x1).*sin(x1);

y1=interp1(x,y,x1);

>> y2=interp1(x,y,x1,'cubic');

>> y3=interp1(x,y,x1,'spline');

>> y4=interp1(x,y,x1,'nearest');

>> plot(x1,[y1;y2;y3;y4],':',x,y,'o',x1,y0)或者plot(x1,[y1',y2',y3',y4'],':',x,y,'o',x1,y0)

2 用不同插值的方法计算sin(x)在pi/2的值

clear

clc

x=0:0.2:pi;%产生包含被插值点的采样点(做成一个向量)

y=sin(x);%求出各采样点对应的样本值

y1=interp1(x,y,pi/2);%用默认的'linear'方法计算sin(pi/2)

y2=interp1(x,y,pi/2,'nearest');%用默认的'linear'方法计算sin(pi/2)

y3=interp1(x,y,pi/2,'cubic');%用三次多项式插值方法计算sin(pi/2)

y4=interp1(x,y,pi/2,'spline');%用三次样条插值方法计算sin(pi/2)

y5=spline(x,y,pi/2);%直接用'spline'方法计算sin(pi/2),功能与y4相同

disp('各种方法的插值结果:')

out=['y1=', num2str(y1)

'y2=',num2str(y2);

'y3=',num2str(y3)]

out2=['y4=',num2str(y4)

'y5=',num2str(y5)]

3求

?????-=-=y x dt dy y

x dt dx 224 的通解

[x,y]=dsolve('Dx=4*x-2*y','Dy=2*x-y', 't')

%P149(3)的求解

[x,y]=dsolve('Dx=a*y','Dy=-b*x ', 'x(0)=x0', 'y(0)=y0', 't')

y=dsolve('Dy=(b*x)/(a*y)', 'y(0)=y0','x');

4

相比较。的数值解，并与精确解求11)(2)0(10,)1(422'++=?????=≤≤+--=t t y y t t t y y

方法一：

%微分方程的M 函数funt.m 文件

function y=funt(t,y)

y=(y^2-t-2)/4/(t+1);

%以下是求解的脚本m 函数，可自由取名，然后在command 窗口调用该函数求解.

ts=[0,10];%自变量的求解区间

y0=2;%初值条件

[t,y]=ode23('funt',ts,y0);%用2,3阶龙格库塔方法求‘funt ’文件里的微分方程

y1=sqrt(t+1)+1;%求精确解对应点上的函数值

[t,y,y1]%以三列的形式显示结果，其中第一列为采样点t 对应的值，第二列为t 对应的数值解，第二列为t 对应的精确解。

方法二：

%或者直接用内建函数编写待求的微分方程，在一个M 脚本文件里执行求解。

funt=inline('(y^2-t-2)/4/(t+1)') %内建函数编写待求的微分方程

ts=[0,10];%自变量的求解区间

y0=2;%初值条件

[t,y]=ode23(funt,ts,y0);%用2,3阶龙格库塔方法求‘funt ’文件里的微分方程

y1=sqrt(t+1)+1;%求精确解对应点上的函数值

[t,y,y1]

5

绘制系统相平面图象。，求解微分方程组，并，，取模型的状态方程表示为

382810)()()()()()()()()()()()(321233223211===-+-=+-=+-=βρσρσσβt x t x t x t x t x

t x t x t x

t x t x t x t x

lorenz

第一种方法建立lorenz 函数模型的状态方程

function xp=lorenz(t,x)%建立lorenz 函数模型的状态方程

xp=[-8/3,0,x(2);0,-10,10;-x(2),28,-1]*x%表明x 是一个三维向量，前面是一个线性方程组的系数矩阵，xp%是一个三维的输出表示x(1),x(2),x(3)的一阶导数向量。

第二种方法建立lorenz 函数模型的状态方程

function xp=lorenz1(t,x)%建立lorenz 函数模型的状态方程

xp=[-8/3*x(1)+x(2)*x(3);-10*x(2)+10*x(3);-x(2)*x(1)+28*x(2)-x(3)];

%求解微分方程

clear;clc;

x0=[0,0,eps]';%三个初值条件构成的向量

[t,x]=ode23('lorenz',[0,100],x0);

[t,x]

plot(t,x),grid,pause%绘制解x(1),x(2),x(3)各自相对于变量t 的图象，按任意键后继续下面的程序 figure;plot3(x(:,1),x(:,2),x(:,3));%绘制解x(1),x(2),x(3)的关于系统的相平面图象

axis([10,40,-20,20,-20,20]);

上机练习题七

班级： 姓名： 学号： 1有一组测量数据如下表所示，数据具有y =x 2+1的变化趋势，用最小二乘法求解y 。

x 1 1.5 2 2.5 3 3.5 4 4.5

5

y -1.4 2.7 3 5.9 8.4 12.2 16.6 18.8 26.2

>> x=[1 1.5 2 2.5 3 3.5 4 4.5 5]'

>> y=[-1.4 2.7 3 5.9 8.4 12.2 16.6 18.8 26.2]'

>> e=[ones(size(x)) x.^2]

>> c=e\y

>> x1=[1:0.1:5]';

>> y1=[ones(size(x1)),x1.^2]*c;

>> plot(x,y,'ro',x1,y1,'k')

2求下列线性方程组的解

（1）????

??????-=????????????????????---129351623114321x x x 解：>> a=[4 1 -1;3 2 -6;1 -5 3];b=[9; -2; 1];x=a\b ％唯一解

x =

2.3830

1.4894

2.0213

（2）????

??????--=????????????????????---121351623334321x x x 解：>> a=[4 -3 3;3 2 -6;1 -5 3];b=[-1;-2;1];x=a\b ％唯一解

x =

-0.4706

-0.2941

（3）????

??????=????????????????-11151231421x x 解：>> a=[4 1;3 2;1 -5];b=[1;1;1];x=a\b ％最小二乘近似解

x =

0.3311

-0.1219

（4）??????????=?????

?????????????????--3211211112111124321x x x x ，求通解 解：>> a=[2 1 -1 1;1 2 1 -1; 1 1 2 1];b=[1;2;3];

>> rank(a),rank([a,b])

ans =

3

ans =

3 ％说明有无穷多解

>> rref([a,b]) ％行最简化

ans =

1.0000 0 0 1.5000 1.0000

0 1.0000 0 -1.5000 0

0 0 1.0000 0.5000 1.0000

％通解为：15.141+-=x x ，425.1x x =，15.043+-=x x

3、求下列矩阵的行列式、逆、特征值和特征向量

（1）????

? ??---351623114

解：>> a=[4 1 -1;3 2 -6;1 -5 3];det(a),inv(a),[v,d]=eig(a)

ans = -94 行列式 ans =

0.2553 -0.0213 0.0426

0.1596 -0.1383 -0.2234 矩阵的逆 0.1809 -0.2234 -0.0532

v =

0.0185 -0.9009 -0.3066

-0.7693 -0.1240 -0.7248 特征向量 -0.6386 -0.4158 0.6170

d =

-3.0527 0 0

0 3.6760 0 特征值

0 0 8.3766

（2）????

? ??--351533134

解：>> a=[4 3 1;3 3 -5;1 -5 3];det(a),inv(a),[v,d]=eig(a)

ans =

-124 行列式

ans =

0.1290 0.1129 0.1452 0.1129 -0.0887 -0.1855 矩阵的逆 0.1452 -0.1855 -0.0242 v = 0.3757 -0.8583 0.3496 -0.6881 -0.0056 0.7255 特征向量 -0.6208 -0.5131 -0.5927

d =

-3.1480 0 0

0 4.6176 0 特征值 0 0 8.5304

(3) ??????

? ??1097591086781075675 解：>> a=[5 7 6 5;7 10 8 7;6 8 10 9;5 7 9 10];det(a),inv(a),[v,d]=eig(a)

ans =

1 行列式

ans =

68.0000 -41.0000 -17.0000 10.0000

-41.0000 25.0000 10.0000 -6.0000 矩阵的逆 -17.0000 10.0000 5.0000 -3.0000

10.0000 -6.0000 -3.0000 2.0000

v =

0.8304 0.0933 0.3963 0.3803

-0.5016 -0.3017 0.6149 0.5286 特征向量 -0.2086 0.7603 -0.2716 0.5520

0.1237 -0.5676 -0.6254 0.5209

d =

0.0102 0 0 0

0 0.8431 0 0

0 0 3.8581 0 特征值 0 0 0 30.2887

(4)

n=5;

tria = 6*eye(n-1); %构造上三角

rown = zeros(1,n-1);

coln = zeros(n,1);