MATLAB程序设计与应用下new.docx
《MATLAB程序设计与应用下new.docx》由会员分享，可在线阅读，更多相关《MATLAB程序设计与应用下new.docx（87页珍藏版）》请在冰豆网上搜索。
MATLAB程序设计与应用下new
MATLAB
程序设计与实例应用
（下）
晋中学院物理与电子工程学院
第六章
MATLAB的数值计算
第六章MATLAB数值计算
6-1多项式的运算
6-1-1多项式的生成和表达
1.多项式的表达
在MATLAB环境下多项式是用向量的形式表达的.向量最右边的元素表示多项式的0阶，向左数依次表示多项式的第1阶、第2阶、第3阶…。
例如多项式
表示为：
[50321]。
2.多项式的生成
语法：
P=ploy（MA）
说明：
1.若MA为方阵，则生成的多项式P为方阵MA的特征多项式。
2.若MA为向量，则向量和多项式满足这样一种关系：
，生成的多项式为：
3.直接输入的方式生成多项式。
例6-1
利用方阵M=[567;891;111213]生成一个多项式（为方阵M的特征多项式）。
程序设计：
>>clear
M=[567;891;111213];
P=poly（M）;%产生多项式的向量表达式
Px=poly2str（P,'x'）;%生成常见的多项式表示形式
P,Px
运行结果：
P=
1.0000-27.000090.000054.0000
Px=
x^3-27x^2+90x+54
例6-2
利用向量A=[2345]生成一个多项式。
程序设计：
>>clear
A=[2345];
P=poly（A）;%产生多项式的向量表达式
Px=poly2str（P,'x'）;%生成常见的多项式表示形式
P,Px
运行结果：
P=
1-1471-154120
Px=
x^4-14x^3+71x^2-154x+120
6-1-2多项式的乘除
语法：
A.c=conv（a,b）
B.[q,r]=decony（c,a）
说明：
1.a、b和c分别是多项式的向量表示形式。
A表示两个多项式的乘积运算，B表示两个多项式的除法运算。
2.q表示除运算的商，r表示除运算的余数。
例6-3
求多项式
和
的乘积M（x）。
程序设计：
>>clear
a=[150];%第一个多项式F（x）
b=[21];%第二个多项式G（x）
c=conv（a,b）;%求两个多项式的乘积
Mx=poly2str（c,'x'）;%用常用的方式表示多项式的积
c,Mx
%end
运行结果：
c=
21150
Mx=
2x^3+11x^2+5x
例6-4
求多项式
和
的除运算D（x）。
程序设计：
>>clear
c=[150];%第一个多项式F（x）
a=[21];%第二个多项式G（x）
[q,r]=deconv（c,a）;%求F（x）/G（x）
Dx=poly2str（q,'x'）;%用常用的方式表示多项式的积
q,r,Dx
%end
运行结果：
q=
0.50002.2500
r=
00-2.2500
Dx=
0.5x+2.25
程序说明：
1.在运行结果中变量q是F（x）除以G（x）的商，而r则是除不尽的余数。
2.运行结果变量Dx表示的商没有加上余数。
6-1-3多项式的求导
语法：
Dp=polyder（p）
说明：
p为向量表示的多项式。
例6-5
求多项式
和
的一阶导数。
我们容易知道以上两个方程的导数手工验算结果为：
F'（x）=2x+5和G'（x）=2
我们看MATLAB的计算结果。
程序设计：
>>clear
f=[150];%第一个多项式F（x）
g=[21];%第二个多项式G（x）
Df=polyder（f）;%求F（x）的导数
Dg=polyder（g）;%求G（x）的导数
Dfx=poly2str（Df,'x'）;
Dgx=poly2str（Dg,'x'）;
Df,Dg,Dfx,Dgx
%end
运行结果：
Df=
25
Dg=
2
Dfx=
2x+5
Dgx=
2
6-1-4多项式的求根
语法：
A.r=roots（p）
说明：
另外还有一种通过先求多项式的伴随矩阵，然后再求特征值的方法也可以求得多项式的根。
例6-6
求多项式
和
的根。
我们容易知道方程的根为：
F（x）为：
x1=0；x2=-5
G（x）为：
x2=-1/2
我们看MATLAB的计算结果。
程序设计：
>>clear
f=[150];%第一个多项式F（x）
g=[21];%第二个多项式G（x）
rf=roots（f）;%求F（x）的根
rg=roots（g）;%求G（x）的根
rf,rg,
%end
运行结果：
rf=
0
-5
rg=
-0.5000
程序说明：
从运行结果我们可以看到F（x）的根为0和-5，G（x）的根为-0.50，这和我们手工验算的结果完全一致。
例6-7
求多项式
和
的根。
程序设计：
>>clear
f=[15-3];%第一个多项式F（x）
g=[2101];%第二个多项式G（x）
rf=roots（f）;%求F（x）的根
rg=roots（g）;%求G（x）的根
rf,rg,
%end
运行结果：
rf=
-5.5414
0.5414
rg=
-1.0000
0.2500+0.6614i
0.2500-0.6614i
程序说明：
1.我们可以看到例6-6和例6-7求得的根和我们手工计算的结果是一致的，其中例6-7多项式的根出现了虚根。
2.我们用求伴随矩阵的方法对例6-7再做一次求根运算：
>>clear
f=[15-3];%第一个多项式F（x）
g=[2101];%第二个多项式G（x）
cf=compan（f）;%求F（x）的伴随矩阵
cg=compan（g）;%求G（x）的伴随矩阵
cf,cg,
%end
运行结果：
cf=
-53
10
cg=
-0.50000-0.5000
1.000000
01.00000
ccf=eig（cf）%求特征根
ccg=eig（cg）%求特征根
运行结果：
ccf=
-5.5414
0.5414
ccg=
-1.0000
0.2500+0.6614i
0.2500-0.6614i
我们看到两种方法的结果是一致的。
6-2数据分析
6-1-1极值、均值、标准差和中位值的计算
语法：
A.Pmax=max（X）
B.Pmin=min（X）
C.Pmean=mean（X）
D.Pstd=std（X）
E.Pmedian=median（X）
说明：
1.max（X）、min（X）、mean（X）、std（X）、median（X）分别是用来求数组或矩阵的最大值、最小值、均值、标准差、中位值。
注意std（X）的定义是
，std（X,1）的定义是
。
2.这里X可以是数组也可以是矩阵。
如果是数组则对整个数组进行计算，如果是矩阵则是分别对矩阵的每个列向量进行运算。
例6-8
对一组随机数组进行均值、方差和中位值的计算。
程序设计：
>>clear
x=randn（1,10）;%生成一个10元素的数组
Pmax=max（x）;%计算数组x的最大值
Pmin=min（x）;%计算数组x的最小值
Pmean=mean（x）;%计算数组x的平均值
Pstd=std（x）;%计算数组x的标准差
Psqu=Pstd^2;%计算数组x的方差
Pmed=median（x）;%计算数组x的中位值
answers=[PmaxPminPmeanPstdPsquPmed];
x,answers
运行结果：
x=
Columns1through6
-0.4326-1.66560.12530.2877-1.14651.1909
Columns7through10
1.1892-0.03760.32730.1746
answers=
1.1909-1.66560.00130.90340.81620.1500
程序说明：
1.注意方差为标准差的平方。
2.需要注意的是MATLAB的变量定义是区分大小写的。
例6-9
对一个随机矩阵进行最大值、最小值、均值、标准差、中位值的计算。
程序设计：
>>clear
X=randn（6,3）;%生成一个6×3的随机矩阵
Pmax=max（X）;%计算矩阵X的最大值
Pmin=min（X）;%计算矩阵X的最小值
Pmean=mean（X）;%计算矩阵X的平均值
Pstd=std（X）;%计算矩阵X的标准差
Pmed=median（X）;%计算矩阵X的中位值
X,Pmax,Pmin,Pmean,Pstd,Pmed;
运行结果：
X=
-0.18671.06680.7143
0.72580.05931.6236
-0.5883-0.0956-0.6918
2.1832-0.83230.8580
-0.13640.29441.2540
0.1139-1.3362-1.5937
Pmax=
2.18321.06681.6236
Pmin=
-0.5883-1.3362-1.5937
Pmean=
0.3519-0.14060.3607
Pstd=
0.99620.84821.2404
6-2-2曲线拟合
语法：
A.P=polyfit（X,Y,N）
B.Yval=polyval（P,X）
说明：
1.polyfit（X,Y,N）根据输入数据X和Y生成一个N阶的拟合多项式。
2.polyval（P,X）根据数据X，用拟合多项式P生成拟合好的数据。
3.这里X和Y构成了一系列数据点的坐标。
例6-10
下列数据为一段时间的日气温平均实际数据，求该数据的拟合方程。
设这组气温数据为：
[181917182022222321212022232322]
程序设计：
>>clear
d=1:
15;%时间
t=[181917182022222321212022232322];%和时间相对应的日平均气温
p=polyfit（d,t,3）;%生成拟合多项式（向量形式），阶次为3
px=poly2str（p,'x'）;%生成拟合多项式（字符形式）
pv=polyval（p,d）;%生成拟合值
p,px
plot（d,t,':
*r',d,pv,'-k','LineWidth',2）
title（'\fontsize{18}\bf数据拟合图像','Color','r'）
xlabel（'\fontsize{14}\rm时间d轴','Color','k'）
ylabel（'\fontsize{14}\rm气温t轴','Color','k'）
gtext（'\fontsize{16}\fontname{宋体}实际曲线','Color','r'）
gtext（'\fontsize{16}\fontname{宋体}拟合曲线','Color','k'）
%end
运行结果：
p=
0.0010-0.05380.964416.4623
px=
0.0010474x^3-0.053823x^2+0.96436x+16.4623
图形如6-1所示：
程序说明：
1.d表示时间第一天、第二天…，t表示每天的气温。
2.本例我们选取的拟合多项式的阶次是3。
3.p为向量形式的拟合多项式，px为字符形式的拟合多项式。
4.plot为绘图函数，它的两个参数分别为图像的横坐标和纵坐标。
4个参数则分别为两条曲线的横坐标和纵坐标。
5.图中的箭头可在绘图窗口中打开insert菜单，选择Arrow然后用鼠标从起始位置拖动到终止位置即可。
例6-11
根据上例中的数据建立的天气温度，拟合方程求5天后的天气气温是多少？
程序设计及运行结果：
>>clear
x=20;
T=0.0010474*x^3-0.053823*x^2+0.96436*x+16.4623
T=
22.5995
程序设计及运行结果：
1.依照上例7-10天的时间序列d，5天后应该是20天，所以令x=20。
2.T即为上例中我们求得的拟合公式，结果求得5天后的日平均气温为22.5995。
6-3数值积分和数值微分
6-3-1微分和积分的数学表达式
对函数f（x）在区间[x1,x2]上求积分在数学上常常表示为
对函数F（x）求导数运算在数学上常常表示为：
6-3-2函数数值积分
命令函数：
1.连续被积函数
quadl（'f',a,b,t）自适应递推牛顿—科西（Newton—Cotes）法求积分
2.离散被积函数
trapz（X,Y）梯形法求数值积分
说明：
1.
1）参数'f'是被积函数表达式字符串或者函数文件名。
2）参数a、b定义函数积分的上下限。
3）参数t定义积分的精度（默认值1.0e-6）。
2.
trapz积分中给出Y相对于X的积分值。
当Y是（m*n）矩阵时，积分对Y的列向量分别进行，得到一个（1*n）矩阵是Y的列向量对应于X的积分结果。
例6-12
设
，用牛顿—科西法求
。
这个积分可以用手工算出来是2，我们看MATLAB的计算结果。
程序设计及运行结果：
>>clear
s=quadl（'sin（x）',0,pi）
s=
2.0000
例6-13
设
，求
。
程序设计及运行结果：
>>clear
y=quadl（'（sin（x）.^2+x.^3）./（x-2.*cos（3.*x））',0,2）
%end
ans=
11.1943
例6-14
用离散数据表示
，积分区间为
分别用离散积分方法对其进行积分。
程序设计及运行结果（如图6-3）：
>>clear
X=0:
0.001:
pi;%离散化积分变量
Y=sin（X）;%生成被积函数的离散序列
Z=trapz（X,Y）;%梯形法求积分
Z
%end
Z=
2.0000
例6-15
设
，求
。
程序设计及运行结果：
>>clear
X=0:
0.001:
2;%离散化积分变量
Y=（sin（X）.^2+X.^3）./（X-2.*cos（3.*X））;
Z=trapz（X,Y）;
Z
%end
Z=
10.9607
程序说明：
输入公式时一定要注意算符都是对数组进行运算的，所以是“.*”“.^”“./“。
6-3-3数值微分
命令函数：
D=diff（Y）
说明：
Y为一组离散序列，D为与Y同维的离散矩阵，则微分运算是对矩阵的每个列向量微分。
例6-16
设函数
和
进行微分运算，并画出计算结果的图像。
程序设计：
>>clear
dx=0.001;
x=-3:
dx:
3;%生成一个从-3到3间隔为0.001的序列
x1=-3:
dx:
3-dx;
y=sin（x）;
z=2.*x.^3+3.*x.^2+x-1;
Y=diff（y）;%对y进行微分运算
DY=diff（y）/dx;%对y进行求导运算
Z=diff（z）;%对z进行微分运算
DZ=diff（z）/dx;%对z进行求导运算
subplot（2,2,1）;%产生两行一列的绘图区间的第一个绘图区间
plot（Y,'-k','LineWidth',1）
title（'\fontsize{12}\fontname{TimeNewRoman}\Deltay（x）','Color','k'）
gridon
subplot（2,2,2）;%产生两行一列的绘图区间的第二个绘图区间
plot（Z,'-k','LineWidth',1）
title（'\fontsize{12}\fontname{TimeNewRoman}\Deltaz（x）','Color','k'）
gridon
subplot（2,2,3）;%产生两行一列的绘图区间的第三个绘图区间
plot（x1,DY,'-k','LineWidth',1）
title（'\fontsize{12}\fontname{TimeNewRoman}y''（x）=cosx','Color','k'）
xlabel（'\fontsize{12}\fontname{TimeNewRoman}x','Color','k'）
ylabel（'\fontsize{12}\fontname{TimeNewRoman}dy/dx','Color','k'）
axis（[-33-11]）
gridon
subplot（2,2,4）;%产生两行一列的绘图区间的第四个绘图区间
plot（x1,DZ,'-k','LineWidth',1）
title（'\fontsize{12}\fontname{TimeNewRoman}z''（x）=6x^2+6x+1','Color','k'）
xlabel（'\fontsize{12}\fontname{TimeNewRoman}x','Color','k'）
ylabel（'\fontsize{12}\fontname{TimeNewRoman}dz/dx','Color','k'）
axis（[-33-2080]）
gridon
%end
运行结果（如图6-2）：
程序说明：
1.输入公式时一定要注意算符都是对数组进行运算的，所以是“.*”“.^”“./”。
2.Y=diff（y）和Z=diff（z）和依此画出的图像，纵坐标是两函数的微分，横坐标是数据序列的值并不是实际图像的横坐标。
由于纵坐标是
，所以纵坐标在数值上是导数的
倍。
3.导数图的纵坐标是
，n的取值范围应该是
-1。
因此
取值是
到
，而微分相应取值是
到
，所以两序列的数目相差1，所以画图时不能用plot（x,diff（y）/dx）而是plot（x1,diff（y）/dx）。
6-3一般非线性方程组的数值解
命令函数：
x=fsolve（'fun',x0,options）
说明：
1.'fun'是我们要求解的方程组，可以直接输入，不过由于方程组比较复杂，一般都先生成一个m文件，然后再进行调用。
2.x0是给出的这个方程组的初值解，我们可以随意地给出。
3.options是命令函数fsolve的参数设置项。
因为fsolve求解的过程本身就是一个优化的过程，而options则是设置优化过程的参数的。
这里对于一般的非线性方程组求解我们常设置为optiomset（'Display','off'），意思为不显示优化过程。
例：
6-17
求方程组
的数值解。
这个方程组的解经过手工演算我们知道解为：
x=4/3，y=-2/3，z=-2/3。
程序设计与运行结果：
1.建立函数文件xyz，并保存：
functionfun=xyz（X）
x=X
（1）;y=X
（2）;z=X（3）;
fun=zeros（3,1）;
fun
（1）=x+y+z;
fun
（2）=x-y+3*z;
fun（3）=3*x+y-z-4;
2.在MATLAB的工作区里求解：
>>X0=[111];
X=fsolve（@xyz,X0,optimset（'Display','off'））
X=
1.3333-0.6667-0.6667
程序说明：
1.在建立一个函数文件时一定要以function作为文件头说明。
2.把建立的函数文件保存到MATLAB默认的文件夹里，如work文件夹。
3.zeros（3,1）用来产生一个3行1列的矩阵，矩阵的所有元素都为0。
4.我们可以看到计算的结果和我们手工计算的结果非常近似，差别是计算精度引起的，fsolve的默认求解精度是0.001。
5.事实上对方程组的求解过程是一个优化过程，optimset（'Display','off'）设置不显示求解过程。
例：
6-18
求非线性方程组
的数值解。
这个方程组手工计算是难以完成的。
程序设计与运行结果：
1.建立函数文件fun，并保存：
functionQ=fun（T）
x=T
（1）;y=T
（2）;z=T（3）;
Q=zeros（3,1）;
Q
（1）=x+y-z;
Q
（2）=cos（x）+y^3+z-5;
Q（3）=4*x+2^y+log（z）-1;
2.在MATLAB的工作区里求解：
>>clear
X=fsolve（@fun,[111],optimset（'Display','off'））
X=
-0.43921.45481.0156
程序说明：
fsolve的默认求解精度是0.001，方程成立的精度为0.001。
6-4微分方程的数值解
6-4-1微分方程的基本形式
的导数等于
，即：
6-4-2一阶常微分方程的求解
语法：
A.[T,Y]=ode45（odefun,tspan,y0）
B.[T,Y]=ode23（odefun,tspan,y0）
说明：
1.ode45（）和ode23（）是两个求解微分方程的命令函数，ode45（）应用最广，在容许误差大的情况下ode23（）的效率要比ode45（）高一些。
2.T和Y分别表示微分方程的解的两个变量的数值序列。
3.odefun是待解微分方程表达式的函数文件名。
4.tspan表示运算的起至时间，是行向量，例如[010]表示运算时间从0开始到10结束。
5.y0初始状态，用列向量表示，其元素分别表示微分方程组各个表达式的初始状态值。
例：
6-19
求解微分方程
并画出解在区间[02*pi]上的图像。
程序设计与运行结果：
1.在程序编辑器中建立待解微分方程表达式的函数文件funx.m，并存入work工作目录：
functionY=funx（x,y）
升级会员 