第六章MATLAB 数值计算Word格式.docx
《第六章MATLAB 数值计算Word格式.docx》由会员分享，可在线阅读，更多相关《第六章MATLAB 数值计算Word格式.docx（26页珍藏版）》请在冰豆网上搜索。
例6-3
求多项式
和
的乘积M（x）。
a=[150];
%第一个多项式F（x）
b=[21];
%第二个多项式G（x）
c=conv（a,b）;
%求两个多项式的乘积
Mx=poly2str（c,'
%用常用的方式表示多项式的积
c,Mx
%end
c=
21150
Mx=
2x^3+11x^2+5x
例6-4
的除运算D（x）。
c=[150];
%第一个多项式F（x）
a=[21];
%第二个多项式G（x）
[q,r]=deconv（c,a）;
%求F（x）/G（x）
Dx=poly2str（q,'
q,r,Dx
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的计算结果。
f=[150];
g=[21];
Df=polyder（f）;
%求F（x）的导数
Dg=polyder（g）;
%求G（x）的导数
Dfx=poly2str（Df,'
Dgx=poly2str（Dg,'
Df,Dg,Dfx,Dgx
Df=
25
Dg=
2
Dfx=
2x+5
Dgx=
6-1-4多项式的求根
A.r=roots（p）
另外还有一种通过先求多项式的伴随矩阵，然后再求特征值的方法也可以求得多项式的根。
例6-6
的根。
我们容易知道方程的根为：
F（x）为：
x1=0；
x2=-5
G（x）为：
x2=-1/2
rf=roots（f）;
%求F（x）的根
rg=roots（g）;
%求G（x）的根
rf,rg,
rf=
0
-5
rg=
-0.5000
从运行结果我们可以看到F（x）的根为0和-5，G（x）的根为-0.50，这和我们手工验算的结果完全一致。
例6-7
f=[15-3];
g=[2101];
%求F（x）的根
%求G（x）的根
-5.5414
0.5414
-1.0000
0.2500+0.6614i
0.2500-0.6614i
1.我们可以看到例6-6和例6-7求得的根和我们手工计算的结果是一致的，其中例6-7多项式的根出现了虚根。
2.我们用求伴随矩阵的方法对例6-7再做一次求根运算：
%第一个多项式F（x）
cf=compan（f）;
%求F（x）的伴随矩阵
cg=compan（g）;
%求G（x）的伴随矩阵
cf,cg,
cf=
-53
10
cg=
-0.50000-0.5000
1.000000
01.00000
ccf=eig（cf）%求特征根
ccg=eig（cg）%求特征根
ccf=
-5.5414
ccg=
-1.0000
我们看到两种方法的结果是一致的。
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
对一组随机数组进行均值、方差和中位值的计算。
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
对一个随机矩阵进行最大值、最小值、均值、标准差、中位值的计算。
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]
d=1:
15;
%时间
t=[181917182022222321212022232322];
%和时间相对应的日平均气温
p=polyfit（d,t,3）;
%生成拟合多项式（向量形式），阶次为3
px=poly2str（p,'
%生成拟合多项式（字符形式）
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轴'
k'
ylabel（'
\fontsize{14}\rm气温t轴'
gtext（'
\fontsize{16}\fontname{宋体}实际曲线'
\fontsize{16}\fontname{宋体}拟合曲线'
p=
0.0010-0.05380.964416.4623
px=
0.x^3-0.x^2+0.96436x+16.4623
图形如6-1所示：
1.d表示时间第一天、第二天…，t表示每天的气温。
2.本例我们选取的拟合多项式的阶次是3。
3.p为向量形式的拟合多项式，px为字符形式的拟合多项式。
4.plot为绘图函数，它的两个参数分别为图像的横坐标和纵坐标。
4个参数则分别为两条曲线的横坐标和纵坐标。
5.图中的箭头可在绘图窗口中打开insert菜单，选择Arrow然后用鼠标从起始位置拖动到终止位置即可。
例6-11
根据上例中的数据建立的天气温度，拟合方程求5天后的天气气温是多少？
程序设计及运行结果：
x=20;
T=0.*x^3-0.*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,trace）自适应递推牛顿—科西（Newton—Cotes）法求积分
2.离散被积函数
trapz（X,Y）梯形法求数值积分
1.
1）参数'
是被积函数表达式字符串或者函数文件名。
2）参数a、b定义函数积分的上下限。
3）参数t定义积分的精度（默认值1.0e-6）；
trace设置是否用图形展示积分过程，1为展示、0不展示。
2.
trapz积分中给出Y相对于X的积分值。
当Y是（m*n）矩阵时，积分对Y的列向量分别进行，得到一个（1*n）矩阵是Y的列向量对应于X的积分结果。
例6-12
设
，用牛顿—科西法求
这个积分可以用手工算出来是2，我们看MATLAB的计算结果。
s=quadl（'
sin（x）'
0,pi）
s=
2.0000
例6-13
，求
y=quadl（'
（sin（x）.^2+x.^3）./（x-2.*cos（3.*x））'
0,2）
ans=
11.1943
例6-14
用离散数据表示
，积分区间为
分别用离散积分方法对其进行积分。
程序设计及运行结果（如图6-3）：
X=0:
0.001:
pi;
%离散化积分变量
Y=sin（X）;
%生成被积函数的离散序列
Z=trapz（X,Y）;
%梯形法求积分
Z
Z=
例6-15
2;
%离散化积分变量
Y=（sin（X）.^2+X.^3）./（X-2.*cos（3.*X））;
Z
10.9607
输入公式时一定要注意算符都是对数组进行运算的，所以是“.*”“.^”“./“。
6-3-3数值微分
D=diff（Y）
Y为一组离散序列，D为与Y同维的离散矩阵，则微分运算是对矩阵的每个列向量微分。
例6-16
设函数
进行微分运算，并画出计算结果的图像。
dx=0.001;
x=-3:
dx:
3;
%生成一个从-3到3间隔为0.001的序列
x1=-3:
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,'
1）
\fontsize{12}\fontname{TimeNewRoman}\Deltay（x）'
gridon
subplot（2,2,2）;
%产生两行一列的绘图区间的第二个绘图区间
plot（Z,'
\fontsize{12}\fontname{TimeNewRoman}\Deltaz（x）'
subplot（2,2,3）;
%产生两行一列的绘图区间的第三个绘图区间
plot（x1,DY,'
1）
\fontsize{12}\fontname{TimeNewRoman}y'
'
（x）=cosx'
xlabel（'
\fontsize{12}\fontname{TimeNewRoman}x'
ylabel（'
\fontsize{12}\fontname{TimeNewRoman}dy/dx'
axis（[-33-11]）
subplot（2,2,4）;
%产生两行一列的绘图区间的第四个绘图区间
plot（x1,DZ,'
\fontsize{12}\fontname{TimeNewRoman}z'
（x）=6x^2+6x+1'
\fontsize{12}\fontname{TimeNewRoman}dz/dx'
axis（[-33-2080]）
运行结果（如图6-2）：
1.输入公式时一定要注意算符都是对数组进行运算的，所以是“.*”“.^”“./”。
2.Y=diff（y）和Z=diff（z）和依此画出的图像，纵坐标是两函数的微分，横坐标是数据序列的值并不是实际图像的横坐标。
由于纵坐标是
，所以纵坐标在数值上是导数的
倍。
3.导数图的纵坐标是
，n的取值范围应该是
因此
取值是
到
，而微分相应取值是
，所以两序列的数目相差1，所以画图时不能用plot（x,diff（y）/dx）。
6-4一般非线性方程组的数值解
x=fsolve（'
fun'
x0,options）
1.'
是我们要求解的方程组，可以直接输入，不过由于方程组比较复杂，一般都先生成一个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（'
））
1.3333-0.6667-0.6667
1.在建立一个函数文件时一定要以function作为文件头说明。
2.把建立的函数文件保存到MATLAB默认的文件夹里，如work文件夹。
3.zeros（3,1）用来产生一个3行1列的矩阵，矩阵的所有元素都为0。
4.我们可以看到计算的结果和我们手工计算的结果非常近似，差别是计算精度引起的，fsolve的默认求解精度是0.001。
5.事实上对方程组的求解过程是一个优化过程，optimset（'
）设置不显示求解过程。
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;
X=fsolve（@fun,[111],optimset（'
-0.43921.45481.0156
fsolve的默认求解精度是0.001，方程成立的精度为0.001。
6-5微分方程的数值解
6-5-1微分方程的基本形式
的导数等于
，即：
6-5-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）
Y=cos（x）;
2.在MATLAB的命令窗口求解，结果如图6-3所示：
[x,y]=ode45（@funx,[02*pi],0）;
升级会员 