例61 求矩阵A的每行及每列的最大和最小元素并求整个矩阵.docx

例61 求矩阵A的每行及每列的最大和最小元素并求整个矩阵.docx

《例61 求矩阵A的每行及每列的最大和最小元素并求整个矩阵.docx》由会员分享，可在线阅读，更多相关《例61 求矩阵A的每行及每列的最大和最小元素并求整个矩阵.docx（20页珍藏版）》请在冰豆网上搜索。

例61求矩阵A的每行及每列的最大和最小元素并求整个矩阵

例6.1求矩阵A的每行及每列的最大和最小元素，并求整个矩阵的最大和最小元素。

A=[13,-56,78;25,63,-235;78,25,563;1,0,-1];

max（A,[],2）%求每行最大元素

min（A,[],2）%求每行最小元素

max（A）%求每列最大元素

min（A）%求每列最小元素

max（max（A））%求整个矩阵的最大元素。

也可使用命令：

max（A（:

））

min（min（A））%求整个矩阵的最小元素。

也可使用命令：

min（A（:

））

例6.2求矩阵A的每行元素的乘积和全部元素的乘积。

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

S=prod（A,2）

prod（S）%求A的全部元素的乘积。

也可以使用命令prod（A（:

））

例6.3求向量X=（1！

，2！

，3！

，…，10！

）。

X=cumprod（1:

10）

例6.4对二维矩阵x，从不同维方向求出其标准方差。

x=[4,5,6;1,4,8]%产生一个二维矩阵x

y1=std（x,0,1）

y2=std（x,1,1）

y3=std（x,0,2）

y4=std（x,1,2）

例6.5生成满足正态分布的10000×5随机矩阵，然后求各列元素的均值和标准方差，再求这5列随机数据的相关系数矩阵。

X=randn（10000,5）;

M=mean（X）

D=std（X）

R=corrcoef（X）

例6.6对下列矩阵做各种排序。

A=[1,-8,5;4,12,6;13,7,-13];

sort（A）%对A的每列按升序排序

-sort（-A,2）%对A的每行按降序排序

[X,I]=sort（A）%对A按列排序，并将每个元素所在行号送矩阵I

例6.7给出概率积分

的数据表如表6.1所示，用不同的插值方法计算f（0.472）。

表6.1概率积分数据表

x

0.46

0.47

0.48

0.49

f（x）

0.4846555

0.4937542

0.5027498

0.5116683

x=0.46:

0.01:

0.49;%给出x，f（x）

f=[0.4846555,0.4937542,0.5027498,0.5116683];

formatlong

interp1（x,f,0.472）%用默认方法，即线性插值方法计算f（x）

interp1（x,f,0.472,'nearest'）%用最近点插值方法计算f（x）

interp1（x,f,0.472,'spline'）%用3次样条插值方法计算f（x）

interp1（x,f,0.472,'cubic'）%用3次多项式插值方法计算f（x）

formatshort

例6.8某检测参数f随时间t的采样结果如表5.1，用数据插值法计算t=2，7，12，17，22，17，32，37，42，47，52，57时的f值。

表5.1检测参数f随时间t的采样结果

t

0

5

10

15

20

25

30

f

3.1025

2.256

879.5

1835.9

2968.8

4136.2

5237.9

t

35

40

45

50

55

60

65

f

6152.7

6725.3

6848.3

6403.5

6824.7

7328.5

7857.6

T=0:

5:

65;

X=2:

5:

57;

F=[3.2015,2.2560,879.5,1835.9,2968.8,4136.2,5237.9,6152.7,...

6725.3,6848.3,6403.5,6824.7,7328.5,7857.6];

F1=interp1（T,F,X）%用线性插值方法插值

F1=interp1（T,F,X,'nearest'）%用最近点插值方法插值

F1=interp1（T,F,X,'spline'）%用3次样条插值方法插值

F1=interp1（T,F,X,'cubic'）%用3次多项式插值方法插值

例6.9设z=x2+y2，对z函数在[0,1]×[0,2]区域内进行插值。

x=0:

0.1:

1;y=0:

0.2:

2;

[X,Y]=meshgrid（x,y）;%产生自变量网格坐标

Z=X.^2+Y.^2;%求对应的函数值

interp2（x,y,Z,0.5,0.5）%在（0.5,0.5）点插值

interp2（x,y,Z,[0.50.6],0.4）%在（0.5,0.4）点和（0.6，0.4）点插值

interp2（x,y,Z,[0.50.6],[0.40.5]）%在（0.5,0.4）点和（0.6，0.5）点插值

%下一命令在（0.5,0.4）,（0.6,0.4）,（0.5,0.5）和（0.6,0.5）各点插值

interp2（x,y,Z,[0.50.6]',[0.40.5]）

例6.10某实验对一根长10米的钢轨进行热源的温度传播测试。

用x表示测量点（米），用h表示测量时间（秒），用T表示测得各点的温度（℃），测量结果如表6.2所示。

表6.2钢轨各点温度测量值

Tx

h

02.557.510

0

9514000

30

884832126

60

6764544841

试用用3次多项式插值求出在一分钟内每隔10秒、钢轨每隔0.5米处的温度。

x=0:

2.5:

10;

h=[0:

30:

60]';

T=[95,14,0,0,0;88,48,32,12,6;67,64,54,48,41];

xi=[0:

0.5:

10];

hi=[0:

10:

60]';

temps=interp2（x,h,T,xi,hi,'cubic'）;

mesh（xi,hi,temps）;

例6.11用一个3次多项式在区间[0，2π]内逼近函数

。

X=linspace（0,2*pi,50）;

Y=sin（X）;

P=polyfit（X,Y,3）%得到3次多项式的系数和误差

X=linspace（0,2*pi,20）;

Y=sin（X）;

Y1=polyval（P,X）

plot（X,Y,':

o',X,Y1,'-*'）

例6.12设

（1）求f（x）+g（x）、f（x）-g（x）。

（2）求f（x）×g（x）、f（x）/g（x）。

f=[3,-5,2,-7,5,6];g=[3,5,-3];g1=[0,0,0,g];

f+g1%求f（x）+g（x）

f-g1%求f（x）-g（x）

conv（f,g）%求f（x）*g（x）

[Q,r]=deconv（f,g）%求f（x）/g（x），商式送Q，余式送r。

例6.13求有理分式的导数。

P=[3,5,0,-8,1,-5];

Q=[10,5,0,0,6,0,0,7,-1,0,-100];

[p,q]=polyder（P,Q）

例6.14已知多项式x4+8x3-10，分别取x=1.2和一个2×3矩阵为自变量计算该多项式的值。

A=[1,8,0,0,-10];%4次多项式系数

x=1.2;%取自变量为一数值

y1=polyval（A,x）

x=[-1,1.2,-1.4;2,-1.8,1.6]%给出一个矩阵x

y2=polyval（A,x）%分别计算矩阵x中各元素为自变量的多项式之值

例6.15仍以多项式x4+8x3-10为例，取一个2×2矩阵为自变量分别用polyval和polyvalm计算该多项式的值。

A=[1,8,0,0,-10];%多项式系数

x=[-1,1.2;2,-1.8]%给出一个矩阵x

y1=polyval（A,x）%计算代数多项式的值

y2=polyvalm（A,x）%计算矩阵多项式的值

例6.16求多项式x4+8x3-10的根。

A=[1,8,0,0,-10];

x=roots（A）

例6.17已知

（1）计算f（x）=0的全部根。

（2）由方程f（x）=0的根构造一个多项式g（x），并与f（x）进行对比。

P=[3,0,4,-5,-7.2,5];

X=roots（P）%求方程f（x）=0的根

G=poly（X）%求多项式g（x）

例6.18设x由[0,2π]间均匀分布的10个点组成，求

的1~3阶差分。

X=linspace（0,2*pi,10）;

Y=sin（X）;

DY=diff（Y）;%计算Y的一阶差分

D2Y=diff（Y,2）;%计算Y的二阶差分，也可用命令diff（DY）计算

D3Y=diff（Y,3）;%计算Y的三阶差分，也可用diff（D2Y）或diff（DY,2）

例6.19设

用不同的方法求函数f（x）的数值导数，并在同一个坐标系中做出f'（x）的图像。

f=inline（'sqrt（x.^3+2*x.^2-x+12）+（x+5）.^（1/6）+5*x+2'）;

g=inline（'（3*x.^2+4*x-1）./sqrt（x.^3+2*x.^2-x+12）/2+1/6./（x+5）.^（5/6）+5'）;

x=-3:

0.01:

3;

p=polyfit（x,f（x）,5）;%用5次多项式p拟合f（x）

dp=polyder（p）;%对拟合多项式p求导数dp

dpx=polyval（dp,x）;%求dp在假设点的函数值

dx=diff（f（[x,3.01]））/0.01;%直接对f（x）求数值导数

gx=g（x）;%求函数f的导函数g在假设点的导数

plot（x,dpx,x,dx,'.',x,gx,'-'）;%作图

例6.20用两种不同的方法求：

先建立一个函数文件ex.m：

functionex=ex（x）

ex=exp（-x.^2）;

然后在MATLAB命令窗口，输入命令：

formatlong

I=quad（'ex',0,1）%注意函数名应加字符引号

I=quadl（'ex',0,1）

例6.21用trapz函数计算：

X=0:

0.01:

1;

Y=exp（-X.^2）;

trapz（X,Y）

例6.22计算二重定积分

（1）建立一个函数文件fxy.m：

functionf=fxy（x,y）

globalki;

ki=ki+1;%ki用于统计被积函数的调用次数

f=exp（-x.^2/2）.*sin（x.^2+y）;

（2）调用dblquad函数求解。

globalki;ki=0;

I=dblquad（'fxy',-2,2,-1,1）

ki

例6.23给定数学函数

x（t）=12sin（2π×10t+π/4）+5cos（2π×40t）

取N=128，试对t从0~1秒采样，用FFT作快速傅立叶变换，绘制相应的振幅-频率图。

N=128;%采样点数

T=1;%采样时间终点

t=linspace（0,T,N）;%给出N个采样时间ti（i=1:

N）

x=12*sin（2*pi*10*t+pi/4）+5*cos（2*pi*40*t）;%求各采样点样本值x

dt=t

（2）-t

（1）;%采样周期

f=1/dt;%采样频率（Hz）

X=fft（x）;%计算x的快速傅立叶变换X

F=X（1:

N/2+1）;%F（k）=X（k）（k=1:

N/2+1）

f=f*（0:

N/2）/N;%使频率轴f从零开始

plot（f,abs（F）,'-*'）%绘制振幅-频率图

xlabel（'Frequency'）;

ylabel（'|F（k）|'）

例6.24用直接解法求解下列线性方程组。

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];

b=[13,-9,6,0]';

x=A\b

例6.25用LU分解求解例6.24中的线性方程组。

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];

b=[13,-9,6,0]';

[L,U]=lu（A）;

x=U\（L\b）

[L,U,P]=lu（A）;

x=U\（L\P*b）

例6.26用QR分解求解例6.24中的线性方程组。

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];

b=[13,-9,6,0]';

[Q,R]=qr（A）;

x=R\（Q\b）

[Q,R,E]=qr（A）;

x=E*（R\（Q\b））

例6.27用Cholesky分解求解例6.24中的线性方程组。

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];

b=[13,-9,6,0]';

R=chol（A）

?

?

?

Errorusing==>chol

Matrixmustbepositivedefinite

Jacobi迭代法的MATLAB函数文件Jacobi.m如下：

function[y,n]=jacobi（A,b,x0,eps）

ifnargin==3

eps=1.0e-6;

elseifnargin<3

error

return

end

D=diag（diag（A））;%求A的对角矩阵

L=-tril（A,-1）;%求A的下三角阵

U=-triu（A,1）;%求A的上三角阵

B=D\（L+U）;

f=D\b;

y=B*x0+f;

n=1;%迭代次数

whilenorm（y-x0）>=eps

x0=y;

y=B*x0+f;

n=n+1;

end

例6.28用Jacobi迭代法求解下列线性方程组。

设迭代初值为0，迭代精度为10-6。

在命令中调用函数文件Jacobi.m，命令如下：

A=[10,-1,0;-1,10,-2;0,-2,10];

b=[9,7,6]';

[x,n]=jacobi（A,b,[0,0,0]',1.0e-6）

Gauss-Serdel迭代法的MATLAB函数文件gauseidel.m如下：

function[y,n]=gauseidel（A,b,x0,eps）

ifnargin==3

eps=1.0e-6;

elseifnargin<3

error

return

end

D=diag（diag（A））;%求A的对角矩阵

L=-tril（A,-1）;%求A的下三角阵

U=-triu（A,1）;%求A的上三角阵

G=（D-L）\U;

f=（D-L）\b;

y=G*x0+f;

n=1;%迭代次数

whilenorm（y-x0）>=eps

x0=y;

y=G*x0+f;

n=n+1;

end

例6.29用Gauss-Serdel迭代法求解下列线性方程组。

设迭代初值为0，迭代精度为10-6。

在命令中调用函数文件gauseidel.m，命令如下：

A=[10,-1,0;-1,10,-2;0,-2,10];

b=[9,7,6]';

[x,n]=gauseidel（A,b,[0,0,0]',1.0e-6）

例6.30分别用Jacobi迭代和Gauss-Serdel迭代法求解下列线性方程组，看是否收敛。

a=[1,2,-2;1,1,1;2,2,1];

b=[9;7;6];

[x,n]=jacobi（a,b,[0;0;0]）

[x,n]=gauseidel（a,b,[0;0;0]）

有了上面这些讨论，下面设计一个求解线性方程组的函数文件line_solution.m。

function[x,y]=line_solution（A,b）

[m,n]=size（A）;

y=[];

ifnorm（b）>0%非齐次方程组

ifrank（A）==rank（[A,b]）

ifrank（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）;%0解

ifrank（A）

disp（'方程组有无穷个解，基础解系为y'）;%非0解

y=null（A,'r'）;

end

end

例6.31求解方程组。

A=[1,-2,3,-1;3,-1,5,-3;2,1,2,-2];

b=[1;2;3];

[x,y]=line_solution（A,b）

例6.32求方程组的通解。

formatrat%指定有理式格式输出

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

b=[1,4,0]';

[x,y]=line_solution（A,b）;

x,y

formatshort%恢复默认的短格式输出

例6.33求

在x0=-5和x0=1作为迭代初值时的零点。

先建立函数文件fz.m：

functionf=fz（x）

f=x-1/x+5;

然后调用fzero函数求根。

：

fzero（'fz',-5）%以-5作为迭代初值

fzero（'fz',1）%以1作为迭代初值

例6.34求下列方程组在（1，1，1）附近的解并对结果进行验证。

首先建立函数文件myfun.m。

functionF=myfun（X）

x=X

（1）;

y=X

（2）;

z=X（3）;

F

（1）=sin（x）+y+z^2*exp（x）;

F

（2）=x+y+z;

F（3）=x*y*z;

在给定的初值x0=1，y0=1，z0=1下，调用fsolve函数求方程的根。

X=fsolve（'myfun',[1,1,1],optimset（'Display','off'））

例6.35求圆和直线的两个交点。

圆：

直线：

先建立方程组函数文件fxyz.m：

functionF=fxyz（X）

x=X

（1）;

y=X

（2）;

z=X（3）;

F

（1）=x^2+y^2+z^2-9;

F

（2）=3*x+5*y+6*z;

F（3）=x-3*y-6*z-1;

再在MATLAB命令窗口，输入命令：

X1=fsolve（'fxyz',[-1,1,-1],optimset（'Display','off'））%求第一个交点

X2=fsolve（'fxyz',[1,-1,1],optimset（'Display','off'））%求第二个交点

例6.36求函数

在区间（-10，1）和（1，10）上的最小值点。

首先建立函数文件fx.m：

functionf=f（x）

f=x-1/x+5;

上述函数文件也可用一个语句函数代替：

f=inline（'x-1/x+5'）

再在MATLAB命令窗口，输入命令：

fminbnd（'fx',-10,-1）%求函数在（-10，-1）内的最小值点和最小值

fminbnd（f,1,10）%求函数在（1，10）内的最小值点。

注意函数名f不用加'

例6.37设

求函数f在（0.5,0.5,0.5）附近的最小值。

建立函数文件fxyz.m：

functionf=fxyz（u）

x=u

（1）;y=u

（2）;z=u（3）;

f=x+y.^2./x/4+z.^2./y+2./z;

在MALAB命令窗口，输入命令：

[U,fmin]=fminsearch（'fxyz',[0.5,0.5,0.5]）%求函数的最小值点和最小值

例6.38求解有约束最优化问题。

首先编写目标函数M文件fop.m。

functionf=fop（x）

f=0.4*x

（2）+x

（1）^2+x

（2）^2-x

（1）*x

（2）+1/30*x

（1）^3;

再设定约束条件，并调用fmincon函数求解此约束最优化问题。

x0=[0.5;0.5];

A=[-1,-0.5;-0.5,-1];

b=[-0.4;-0.5];

lb=[0;0];

option=optimset;option.LargeScale='off';option.Display='off';

[x,f]=fmincon（'fop',x0,A,b,[],[],lb,[],[],option）

例6.39设有初值问题：

y（0）=2

试求其数值解，并与精确解相比较（精确解为y（t）=

）。

（1）建立函数文件funt.m。

functiony=funt（t,y）

y=（y^2-t-2）/4/（t+1）;

（2）求解微分方程。

t0=0;tf=10;

y0=2;

[t,y]=ode23（'funt',[t0,tf],y0）;%求数值解

y1=sqrt（t+1）+1;%求精确解

plot（t,y,'b.',t,y1,'r-'）;%通过图形来比较

例6.40已知一个二阶线性系统的微分方程为：

其中a=2，绘制系统的时间响应曲线和相平面图。

建立一个函数文件sys.m：

functionxdot=sys（t,x）

xdot=[-2*x

（2）;x

（1）];

取t0=0，tf=20，求微分方程的解：

t0=0;tf=20;

[t,x]=ode45（'sys',[t0,tf],[1,0]）;

[t,x]

subplot（1,2,1）;plot（t,x（:

2））;%解的曲线，即t-x

subplot（1,2,2）;plot（x（:

2）,x（:

1））%相平面曲线，即x-x’

axisequal

例6.41设

将X转化为稀疏存储方式。

X=[2,0,0,0,0;0,0,0,0,0;0,0,0,5,0;0,1,0,0,-1;0,0,0,0,-5];

A=sparse（X）

例6.42根据表示稀疏矩阵的矩阵A，产生一个稀疏存储方式矩阵B。

A=[2,2,1;3,1,-1;4,3,3;5,3,8;6,6,12];

B=spconvert（A）

例6.43求下列三对角线性代数方程组的解。

B=[1,2,0;1,4,3;2,6,1;1,6,4;0,1,2];%产生非0对角元素矩阵

d=[-1;0;1];%产生非0对角元素位置向量

A=spdiags（B,d,5,5）%产生稀疏存储的系数矩阵

f=[0;3;2;1;5];%方程右边参数向量

x=（inv（A）*f）'%求解