修订版第章MATLAB数值计算例题源程序汇总精Word格式.docx
《修订版第章MATLAB数值计算例题源程序汇总精Word格式.docx》由会员分享，可在线阅读，更多相关《修订版第章MATLAB数值计算例题源程序汇总精Word格式.docx（19页珍藏版）》请在冰点文库上搜索。
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）
spline'
）%用3次样条插值方法计算f（x）
cubic'
）%用3次多项式插值方法计算f（x）
formatshort
例6.8某检测参数f随时间t的采样结果如表6.2，用数据插值法计算t=2，7，12，17，22，17，32，37，42，47，52，57时的f值。
表6.2检测参数f随时间t的采样结果
t
5
10
15
20
25
30
f
3.1025
2.256
879.5
1835.9
2968.8
4136.2
5237.9
35
40
45
50
55
60
65
6152.7
6725.3
6848.3
6403.5
6824.7
7328.5
7857.6
T=0:
5:
65;
X=2:
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,'
）%用最近点插值方法插值
）%用3次样条插值方法插值
）%用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.3钢轨各点温度测量值
Tx
h
02.557.510
9514000
884832126
6764544841
试用3次多项式插值求出在一分钟内每隔10秒、钢轨每隔0.5米处的温度。
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:
temps=interp2（x,h,T,xi,hi,'
）;
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）;
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计算该多项式的值。
%多项式系数
x=[-1,1.2;
2,-1.8]%给出一个矩阵x
y1=polyval（A,x）%计算代数多项式的值
y2=polyvalm（A,x）%计算矩阵多项式的值
例6.16求多项式x4+8x3-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）;
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:
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命令窗口，输入命令：
I=quad（'
ex'
0,1）%注意函数名应加字符引号
I=quadl（'
0,1）
例6.21用trapz函数计算：
X=0:
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函数求解。
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中的线性方程组。
[L,U]=lu（A）;
x=U\（L\b）
[L,U,P]=lu（A）;
x=U\（L\P*b）
例6.26用QR分解求解例6.24中的线性方程组。
[Q,R]=qr（A）;
x=R\（Q\b）
[Q,R,E]=qr（A）;
x=E*（R\（Q\b））
例6.27用Cholesky分解求解例6.24中的线性方程组。
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）
G=（D-L）\U;
f=（D-L）\b;
y=G*x0+f;
y=G*x0+f;
例6.29用Gauss-Serdel迭代法求解下列线性方程组。
在命令中调用函数文件gauseidel.m，命令如下：
[x,n]=gauseidel（A,b,[0,0,0]'
例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;
有了上面这些讨论，下面设计一个求解线性方程组的函数文件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%方程组有无穷多个解，基础解系
原方程组有无穷个解，特解为x，其齐次方程组的基础解系为y'
y=null（A,'
r'
else
方程组无解'
%方程组无解
x=[];
else%齐次方程组
原方程组有零解x'
x=zeros（n,1）;
%0解
ifrank（A）<
n
方程组有无穷个解，基础解系为y'
%非0解
例6.31求解方程组。
A=[1,-2,3,-1;
3,-1,5,-3;
2,1,2,-2];
b=[1;
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作为迭代初值
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）
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（'
））%求第一个交点
X2=fsolve（'
[1,-1,1],optimset（'
））%求第二个交点
例6.36求函数
在区间（-10,1）和（1,10）上的最小值点。
首先建立函数文件fx.m：
functionf=f（x）
上述函数文件也可用一个语句函数代替：
x-1/x+5'
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（'
[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='
option.Display='
[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，求微分方程的解：
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;
1];
%产生非0对角元素位置向量
A=spdiags（B,d,5,5）%产生稀疏存储的系数矩阵
f=[0;
5];
%方程右边参数向量
x=（inv（A）*f）'
%求解
升级会员 