例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）'%求解
升级会员 