第四节数值计算

第四节数值计算第一页，共五十三页，编辑于2023年，星期一学习目标掌握数据分析函数库（datafun）中常用函数的格式和功能掌握矩阵操作函数库（elmat）中常用函数的格式和功能掌握多项式函数库（polyfun）中常用函数的格式和功能能够利用matlab帮助文件来查询函数的格式和使用第二页，共五十三页，编辑于2023年，星期一本节内容掌握统计分析中常用函数的应用掌握复数运算中常用函数的格式和使用掌握多项式运算和矩阵运算掌握曲线拟合和数据插值的应用了解matlab在微积分方面的应用第三页，共五十三页，编辑于2023年，星期一§4.1

复数运算1.复数的表示形式其中：c是复数，a为实部，b为虚部，i是虚数单位其中：z是矢量长度，是弧度另一种表示形式：第四页，共五十三页，编辑于2023年，星期一两种表示形式的换算关系实轴虚轴第五页，共五十三页，编辑于2023年，星期一常用复数计算函数conj(x)：求复数x的共轭数real：求复数的实部imag：求复数的虚部abs：求复数的模值angle：求复数的相位角isreal：判断变量是否为实数，是则返回1第六页，共五十三页，编辑于2023年，星期一§4.2数据统计分析求最大值max和最小值min可查阅datafun函数库中相应的函数。>>docdatafun调用格式如下：

y=max(x)：返回向量x的最大值并存入y；若x中包含复数元素，则按模取最大值[y,I]=max(x)：返回向量x的最大值并存入y；最大值的序号存入I第七页，共五十三页，编辑于2023年，星期一

Y=max(A)：返回矩阵A的每一列上的最大值[Y,U]=max(A)：返回值Y和U均为行向量，Y表示A的每列的最大值，U表示每列最大值的行号

Y=max(A,[],dim)：dim取1时，表示求每一列上的最大值；dim取2时，表示求矩阵A每一行上的最大值

求最小值的函数是min，其用法和max完全相同第八页，共五十三页，编辑于2023年，星期一例：求3阶魔方矩阵中各列和各行元素的最大值，

并求整个矩阵的最大值和最小值。程序如下：

x=magic(3)

max(x)%求矩阵x中各列元素中的最大值

max(x,[],2)%求x中各行元素中的最大值

max(max(x))%求x的最大值9

min(min(x))%求矩阵x的最小值1x=816357492ans=897ans=879第九页，共五十三页，编辑于2023年，星期一数据统计常用函数max求矩阵的最大元素min求最小元素mean求平均值median求中值sort各列元素按升序排列prod各元素的积sortrows各行元素按升排列std求矩阵的标准方差sum求矩阵各列元素的和第十页，共五十三页，编辑于2023年，星期一例:分别求该数组x中各列和各行元素的和。

>>x=[456;148]

x=456

148>>y=sum(x,1)%求各列的和y=594>>y=sum(x,2)%求各行的和y=1513注意：

Y=sum(X,DIM)

％DIM＝1时，按列操作；

DIM＝2时，按行操作第十一页，共五十三页，编辑于2023年，星期一§4.3多项式运算注：行向量中必须包括具有零系数的项，并按x的幂降序排列，若多项式中缺某幂次项，则认为该项的系数为0。

多项式的表示方法是将n阶多项式存储在长度为n+1的行向量p中，行向量p的元素为多项式的系数，并按x的幂的降序排列。第十二页，共五十三页，编辑于2023年，星期一表示方法直接输入多项式的系数向量利用poly指令P=poly(x)功能：以x为具有n个元素的列向量，作为多项式的根来构造多项式，并将多项式的系数赋给向量p第十三页，共五十三页，编辑于2023年，星期一例：把多项式a(x)=x3+9x2+5x-8与b(x)=x3-3x2+28

相加，而后与c(x)=(x-1)(x-2)(x-3)相减。程序如下：

a=[195-8];b=[1-3028];c=poly([123])d1=a+bd2=a-cds2=poly2str(d2,’s’)结果如下：c=1-611-6d1=26520即：c(x)=(x-1)(x-2)(x-3)=x3-6x2+11x-6d2=015-6-2ds2=

15s^2-6s-2第十四页，共五十三页，编辑于2023年，星期一多项式求根令p(x)=0,多项式转换为x的一元n次方程，该方程有n个根，这些根可能为实根也可能是若干对共轭复根。x=roots(P)功能：返回多项式P的全部根，并以列向量形式赋给x.按惯例规定，多项式是行向量，根是列向量。第十五页，共五十三页，编辑于2023年，星期一例：已知f(x)=3x5+4x3-5x2-7.2x+5。

1)计算f(x)=0的全部根；

2）由方程f(x)=0的根构造一个多项式g(x).程序如下：

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

%求方程f(x)=0的根

x=roots(P);%构造多项式g(x)

g=poly(x);结果如下：g=

1.00000.00001.3333-1.6667-2.40001.6667x=-0.3046+1.6217i-0.3046-1.6217i-1.00661.01900.5967第十六页，共五十三页，编辑于2023年，星期一多项式的运算函数c=conv(a,b)多项式a和b的乘积

[q,r]=deconv(b,a)多项式b除以a，q为商，r为余数

[r,p,k]=residue(b,a)有理分式展开思考：residue和residuez之间有什么区别？第十七页，共五十三页，编辑于2023年，星期一程序如下：

>>a=[1234],b=[14916]>>c=conv(a,b)>>[q,r]=deconv(c,a)例：求两多项式相乘结果如下：

c=[162050758464]q＝[14916]r=[0000]第十八页，共五十三页，编辑于2023年，星期一例：将下列有理多项式转化为部分分式.程序如下：>>b=[-25],a=[163];>>[r,p,k]=residue(b,a)即：结果如下：r=-3.2454p=-5.4495k=[]1.2454-0.5505第十九页，共五十三页，编辑于2023年，星期一多项式求导数：polyderr=polyder(P)求多项式P的导数r=polyder(P,Q)求多项式P和Q乘积的导数[r,s]=polyder(P,Q)

求多项式P/Q的导数，导函数的分子存入

r，分母存入s第二十页，共五十三页，编辑于2023年，星期一程序如下：

>>p=[43-21],polyder(P)例：求f(x)=4x3+3x2-2x+1的导数。结果如下：

ans=126-2即：df(x)/x=12x2+6x-2例：求有理分式1/(x2+5)的导数。程序如下：

>>p=1,q=[105],[r,s]=polyder(p,q)结果如下：

r=-20s=1010025即：第二十一页，共五十三页，编辑于2023年，星期一多项式估值：polyvaly=polyval(p,x)求多项式P在点x处的值若x为向量或矩阵，则对其中的每个元素求多项式的值，返回值为与自变量同型Y=polyvalm(P,A)求矩阵多项式的值

以方阵A为自变量求多项式P的值如：p(x)=x2+2x+3,则polyvalm(P,A)表示x*x+2*x+3*eye(size(A))第二十二页，共五十三页，编辑于2023年，星期一程序如下：

>>p=poly([1234]),polyval(p,8)例：求f(x)=(x-1)(x-2)(x-3)(x-4)在x=8处的值。例：当x=[12;34]时，求p(x)=x2+2x+3的值。程序如下：

>>clear;>>P=[123],x=[12;34]>>pa=x.^2+2*x+3;paa=polyval(p,x)>>pm=x^2+2*x+3>>pmm=polyvalm(p,x)第二十三页，共五十三页，编辑于2023年，星期一

多项式拟合指已知变量x、y之间的函数关系为现通过测量数据求多项式系数§4.4曲线拟合

通过试验获得一组测量数据，则求出变量之间的函数关系的过程就叫曲线拟合。曲线拟合的方法有多种如：线性拟合、多项式拟合、指数拟合等。第二十四页，共五十三页，编辑于2023年，星期一

p=polyfit(x,y,n)求数组x,y所给数据的n阶拟合多项式系数向量p1）先用polyfit函数进行多项式拟合，产生一个指定次数的最优拟合多项式系数向量；2）再用polyval求多项式拟合的取值。

为保证较好的拟合效果，多项式系数要取得适当。过低，可能残差比较大；过高，拟合模型将包含噪声影响，通常保证n<m。观察数据图形曲线，有助于适当阶数的确定。第二十五页，共五十三页，编辑于2023年，星期一例：已知实验数据，求4次拟合多项式系数，并计算f(0.4)的值。x=[12345]y=[5.5,43.1,128,290.7,498.4]p=polyfit(x,y,4)

%采用四阶多项式拟合曲线x1=1:0.1:5;y1=polyval(p,x1)plot(x,y,’*b’,x1,y1,’r’)拟合曲线为：-2.6375s^4+31.4583s^3-99.1625s^2+154.4417s-78.6第二十六页，共五十三页，编辑于2023年，星期一数据插值

interp1：一维插值

interp2：二维插值

interp3：三维插值

interpft：一维快速傅里叶插值方法

spline：立体样条插值插值就是在已知数之间寻找估计值的过程。第二十七页，共五十三页，编辑于2023年，星期一

x,y是数据的向量，其长度必须相同；xi是要插值的数据点向量；method为插值方法，可为nearest（最近点插值）、cubic（分段三次Hermite插值）、linear（线性）、spline（三次样本插值）。

interp1(x,y,xi,method)第二十八页，共五十三页，编辑于2023年，星期一例：t表示1900～1990年的人口普查，p表示当年的人口数量，现估计1975年的人口数量。clf%清除图形窗口中的内容t=1900:10:1999;p=[75.99591.972105.711123.203131.669...150.697179.323203.212226.505249.633]interp1(t,p,1975)x=1900:1:2000;%为更详细地估计人口，可插入y=interp1(t,p,x,‘spline’);

%1900-2000年的人口数量plot(t,p,‘o’,x,y)

%画图直观地分析人口增长与时间的关系第二十九页，共五十三页，编辑于2023年，星期一第三十页，共五十三页，编辑于2023年，星期一第三十一页，共五十三页，编辑于2023年，星期一第三十二页，共五十三页，编辑于2023年，星期一§4.5矩阵运算det(A)A的行列式eig(A)A的特征值norm(A,1)A的1范数norm(A)A的2范数norm(A,inf)A的无穷范数norm(A,‘fro’)A的F范数rank(A)A的秩inv(A)A的逆矩阵trace(A)A的迹：矩阵对角线元素之和第三十三页，共五十三页，编辑于2023年，星期一线性方程组的解法它可以表示为,其中线性方程组第三十四页，共五十三页，编辑于2023年，星期一利用左除运算符leftdivisionoperation

x=A\B若方程是xC=d形式，就用“右除”，即x=d/C。采用逆阵法x=inv(A)*B但该法只适用于当A非奇异、维数较低且条件较好时，使用范围有限。上述方程组的解为：第三十五页，共五十三页，编辑于2023年，星期一例：分别采用逆阵法和左除法求解。

已知，求，使解1：逆阵法A=[101;210;-32-5];b=[12-1]’;X=inv(A)*b

解2：直接法ticX=A\btoc所以，左除法速度快，而且精度高。Elapsedtimeis36.672815seconds.Elapsedtimeis3.346355seconds.第三十六页，共五十三页，编辑于2023年，星期一§4.6数据计算求函数极值

一元函数在区间(x1,x2)中的极小值

求多元函数的极小值，指令fminsearch

求极大值：求-f(x)的极小值图形法：指令[xx,yy]=ginput(1)

画出函数图形，观察其极值点，并在极值点附近进行局部放大，在指令窗运行ginput指令，鼠标变为十字型，对准极值点点击即可得到其数组。[x,fval,exitflag,output]=fminbnd(fun,x1,x2,options)第三十七页，共五十三页，编辑于2023年，星期一例如：求函数[0,3]之间的极大值。程序如下：

f=‘exp(-x).*sin(x)’;%字符串类型fx=strcat(‘-’,f);%字符串连接符号xmin=fminbnd(fx,0,3)%极小值x=xmin;ymin=eval(f)%字符串演算指令结果：极大值为（0.7854,0.3224)第三十八页，共五十三页，编辑于2023年，星期一代数方程的解：即f(x)=0的解

求零点fzero解非线性方程组：指令fsolvefzero是根据函数是否穿越横轴来决定零点的。因此它无法确定函数曲线仅触及横轴而不穿越的那些零点。如：(x-1)2sinx中x=1的零点；[t,Y]=ode45(odefun,tspan,y0)

采用四阶龙格－库塔法数值积分法解常数微分方程的解：指令ode45第三十九页，共五十三页，编辑于2023年，星期一例如：求方程组在（0.5,0.5)附近的数值解。解：1）建立函数文件myfun.mfunctionq=myfun(p)x=p(1)；y=p(2);q(1)=x-0.6*sin(x)-0.3*cos(y);q(2)=y-0.6*cos(x)+0.3*sin(y);2）在给定的初值下，调用fsolve函数x=fsolve('myfun',[0.5,0.5])q=myfun(x)第四十页，共五十三页，编辑于2023年，星期一§4.7

数值微积分

在MATLAB数值计算中，由于数值精度有限，故不能表示无穷小量，也不能准确描述一个数的邻域。所以：

不要企图借助数值计算求取极限求微积分最好采用“符号计算”利用数值计算求积分只能得到近似值第四十一页，共五十三页，编辑于2023年，星期一基本思想：将整个积分区间[a，b]分成n个子区间[xi,xi+1]，i=1,2,…,n，其中x1=a，xn+1=b。这样求定积分问题就分解为求和问题。数值积分Sx=sum(X)沿列方向求和

St=trapz(x,y)采用梯形法沿列方向求函数y关于自变量x的积分第四十二页，共五十三页，编辑于2023年，星期一例：求积分。程序如下：x=0:pi/100:pi;y=sin(x);

s=trapz(x,y);

％通用格式，用于“非等间隔采样”

或s=pi/100*trapz(y)

％表示由x,y所绘出的折线下的面积

利用符号计算方法

symsx

s=int(sin(x),x,0,pi)显示结果为s=-2显示结果s=-1.9998第四十三页，共五十三页，编辑于2023年，星期一S1＝quad(fun,a,b,tol)

采用递推自适应Simpson法计算一重积分S1=quadl(fun,a,b,tol)

采用递推自适应Lobatto法求数值的一重积分S2＝blquad(fun,xmin,xmax,ymin,ymax,tol)

二重数值积分三重数值积分函数：triplequad注：1）fun是被积函数，可为字符串、内联对象、匿名函数和M函数文件的函数句柄；2）tol用来控制绝对误差，默认值10-6；3）内向外进行积分的上下限分别是由xmin，xmax,

ymin,ymax传递。第四十四页，共五十三页，编辑于2023年，星期一例：已知螺旋线

求它的弧长

functionf=len(t)f=sqrt(cos(t).^2+4*sin(2*t).^2+1);解：1)建立一个求弧长的被积函数的M-文件len.m2)进入MATLAB环境，在命令窗键入

>>t=0:pi/100:2*pi;>>plot3(sin(t),cos(t),t)

%画出函数的图形

>>length=quad('len',0,2*pi)第四十五页，共五十三页，编辑于2023年，星期一显示结果：

length=11.4609length=quad(@(t)(sqrt(cos(t).^2+4*sin(2*t).^2+1)),0,2*pi)或者：length=quad('(sqrt(cos(t).^2+4*sin(2*t).^2+1))',0,2*pi)另一种求解方法：第四十六页，共五十三页，编辑于2023年，星期一例：求解：f='y*sin(x)+x*cos(y)‘;

或者用f=@(x,y)(y*sin(x)+x*cos(y))

s=dblquad(f,pi,2*pi,0,pi)另解：symsxyf=y*sin(x)+x*cos(y)s=vpa(int(int(f,x,pi,2*pi),y,0,pi))

显示结果为：s=-9.8696显示结果为s=-9.8696044010893586188344909998761第四十七页，共五十三页，编辑于2023年，星期一

所谓的符号计算是指：解算数学表达式方程时，无须事先对变量赋值，即不是在离散化的数值点上进行，而是凭借一系列恒等式和数学定理，通过推理和演绎，获得解析结果，并以标准的符号形式来表示。这种计算建立在数值完全准确表达和推演严格解析的基础之上，因此所得结果是完全准确的。§4.8

符号计算第四十八页，共五十三页，编辑于2023年，星期一定义：利用函数sym或syms

符号变量名＝sym（‘符号字符串’

）

syms符号变量名1变量名2…变量名n注意：sym每次只能定义一个变量，而syms可同时定义多个变量，且各变量名之间只能用空格分隔；例：>>M=sym('[a,b;c,d]')或者采用：>>symsabcd>>M=[a,b;c,d]第四十九页，共五十三页，编辑于2023年，星期一注明：判断变量类型用指令class，其中