数值计算功能

关于数值计算功能第1页，讲稿共81页，2023年5月2日，星期三数值运算的功能创建矩阵矩阵运算多项式运算线性方程组数值统计线性插值函数优化微分方程的数值解第2页，讲稿共81页，2023年5月2日，星期三一、Matlab的基本计算功能函数名称函数功能函数名称函数功能abs(x)取绝对值sign(x)符号函数angle(z)复数的相角rem(x,y)x/y取余sqrt(x)开平方gcd(x,y)最大公因数real(z)复数的实部lcm(x,y)最小公倍数imag(z)复数的虚部exp(x)自然指数conj(z)共轭复数pow2(x)2的指数1、常用基本数学函数P38第3页，讲稿共81页，2023年5月2日，星期三以2为底的对数log2(x)朝零方向取整fix(x)自然对数log(x)四舍五入取整round(x)以10为底的对数log10(x)将实数化为分数rat(x)求余值为x-y.*floor(x./y)mod(x,y)朝正无穷大方向取整ceil(x)求向量长度length(x)朝负无穷方向取整floor(x)函数功能函数名称函数功能函数名称2、Matlab常用的三角函数P38有：sin(x),cos(x),tan(x),sind(x),tand(x)asin(x),acos(x),atan(x)等第4页，讲稿共81页，2023年5月2日，星期三二、命令行的基本操作创建矩阵的方法(1)直接输入法规则：矩阵元素必须用[]括住矩阵元素必须用逗号或空格分隔在[]内矩阵的行与行之间必须用分号分隔

第5页，讲稿共81页，2023年5月2日，星期三

矩阵元素可以是任何matlab表达式，可以是实数，也可以是复数，复数可用特殊函数i，j输入

a=[1,2,3;4,5,6]

矩阵元素a=123456x=[2pi/2;sqrt(3)3+5i]第6页，讲稿共81页，2023年5月2日，星期三符号的作用逗号和分号的作用

逗号和分号可作为指令间的分隔符，matlab允许多条语句在同一行出现。

分号如果出现在指令后，屏幕上将不显示结果，但存储在工作空间中。当一个指令或矩阵太长时，可用•••续行第7页，讲稿共81页，2023年5月2日，星期三用于选出矩阵指定行、列及元素。a=A(:,2:4)或a=A(:,[2,3,4])冒号的作用：用于生成等间隔的向量，默认间隔为1。i=3:2:7循环语句forn=1:10第8页，讲稿共81页，2023年5月2日，星期三(2)用matlab函数创建矩阵(P19-20)空阵[]—matlab允许输入空阵，当一项操作无结果时，返回空阵。rand(m,n)——随机矩阵或rand(n)eye(m,n)——单位矩阵或eye(n)zeros(m,n)——全部元素都为0的矩阵ones(m,n)——全部元素都为1的矩阵linspace(a,b,n))—生成a~b之间n个数值线性分布的向量logspace(a,b,n))—生成10a~10b之间按对数等分的n个元素的向量第9页，讲稿共81页，2023年5月2日，星期三r=eye(3)

r=100010001r=eye(3,4)r=100001000010r=eye(4,3)r=100010001000第10页，讲稿共81页，2023年5月2日，星期三

还有伴随矩阵、稀疏矩阵、魔方矩阵、对角矩阵、范德蒙等矩阵的创建，就不一一介绍了。注意：matlab严格区分大小写字母，因此a与A是两个不同的变量。

第11页，讲稿共81页，2023年5月2日，星期三(3)用M文件创建矩阵适用于较大的矩阵且需经常调用的矩阵。%mymatrix.mcreationofmatrixJZJZ=[1,2,3,4,5;6,7,8,9,10;11,12,13,14,15];取mymatrix名保存若JZ已由其它文件运行后生成而需经常调用，可用save指令保存。savemymatrixJZ(4)从外部数据文件调入矩阵loadmymatrixJZ第12页，讲稿共81页，2023年5月2日，星期三三、矩阵的修改与操作直接修改可用键找到所要修改的矩阵，用键移动到要修改的矩阵元素上即可修改。指令修改可以用A(,)=

来修改。(1)矩阵的修改(5)多维矩阵第13页，讲稿共81页，2023年5月2日，星期三例如a=[120;305;789]a=120305789a(3,3)=0a=120305780还可以用函数subs、find函数修改。第14页，讲稿共81页，2023年5月2日，星期三(2)矩阵的子阵

矩阵的子阵可以通过标量、向量、冒号标志来引用和赋值A(v,w),v,w中任何一个可以是冒号“：”、标量、向量常见A(v,w)形式有：A(i,j)、A(:,j)、A(:,j:k)、A(i:j,k:h)、A([i,j],[k,h])、

A(i:j)等第15页，讲稿共81页，2023年5月2日，星期三(3)矩阵的操作空阵可利用[]清除矩阵中部分行或列来改变维数,其作用与借助向量标识得到的矩阵子块相同。例：已知A=[1,2,3,4;5,6,7,8;9,10,11,12]，要消去A中的第1列与第4列得到B阵。B=A;B(:,[1,4])=[],与B=A(:,[2,3])或B=A(:,2:3)结果相同

23671011第16页，讲稿共81页，2023年5月2日，星期三矩阵的扩展a.利用矩阵标识块赋值指令扩展X(m1:m2,n1:n2)=A生成新阵X，矩阵除了赋值阵A和已存在的元素外，其余为0b.利用方括号和小矩阵生成大矩阵例：利用A=[1,2,7,8;-1,19,7,10]生成4×4的B阵，使其第1、2行，第1至3列元素取A阵第2列至第4列元素，第3、4行元素全为1，其余为0。第17页，讲稿共81页，2023年5月2日，星期三法一：B=zeros(4,4);B(1:2,1:3)=A(:,2:4);B(3:4,:)=ones(2,4)法二：B=[A(:,2:4)，zeros(2,1);ones(2,4)]

278019710011111111或B(1:2,1:3)=A(:,2:4);B(3:4,1:4)=ones(2,4)第18页，讲稿共81页，2023年5月2日，星期三矩阵的结构变换P91-94rot90(A,k):逆时针旋转fliplr:左右翻flipud:上下翻diag:抽取主对角线或生成对角阵tril(A,k):抽取主下三角(k=0,+1,-1,…)triu(A,k):抽取主上三角(k=0,+1,-1,…)B=reshape(A,m,n）根据A阵重组为m×n的B阵第19页，讲稿共81页，2023年5月2日，星期三矩阵加、减（＋,－）运算规则：与线性代数运算规则相同四、矩阵运算2.矩阵乘（）运算规则：与线性代数运算规则相同第20页，讲稿共81页，2023年5月2日，星期三A=[1,3,pi,i;6,8,3+i,5]；B=[1,2,3,4;5,6,7,8]；D=A+B，C=3+BD=2.00005.00006.14164.0000+1.0000i11.000014.000010.0000+1.0000i13.0000a=[123;456;780];b=[1;2;3];c=a*bc=143223第21页，讲稿共81页，2023年5月2日，星期三

a^p——a自乘p次幂方阵>1的整数3.矩阵乘方——a^n,a^p,p^a对于p的其它值,计算将涉及特征值和特征向量，如果p是矩阵，a是标量a^p使用特征值和特征向量自乘到p次幂；如a,p都是矩阵，a^p则无意义。第22页，讲稿共81页，2023年5月2日，星期三inv——矩阵求逆rank——求矩阵的秩det——行列式的值eig——矩阵的特征值diag——对角矩阵’——矩阵转置sqrt——矩阵开方

4.矩阵的其它运算

(P96)第23页，讲稿共81页，2023年5月2日，星期三

数组运算指元素对元素的算术运算，与通常意义上的由符号表示的线性代数矩阵运算不同数组加减(.+,.-)

a.+ba.-b5.矩阵的数组运算

对应元素相加减（与矩阵加减等效）第24页，讲稿共81页，2023年5月2日，星期三2.数组乘除(，./，.\）ab——a，b两数组必须有相同的行和列两数组相应元素相乘。a=[123;456;789];b=[246;135;7910];a.*bans=281841530497290

第25页，讲稿共81页，2023年5月2日，星期三a./b=b.\a—a的元素被b的对应元素除a.\b=b./a—b的元素被a的对应元素除例:

a=[123];b=[456];c1=a.\b;c2=b./ac1=4.00002.50002.0000c2=4.00002.50002.0000第26页，讲稿共81页，2023年5月2日，星期三3.数组乘方(.^)—元素对元素的幂例:a=[123];b=[456];z=a.^2z=1.004.009.00z=a.^bz=1.0032.00729.00第27页，讲稿共81页，2023年5月2日，星期三五、关系运算、逻辑运算及其函数

关系运算和逻辑运算均按照数组运算的规则和定义进行的。1、关系运算（2）运算规则>,<,>=,>=,==,~=（1）关系运算符

a.两个标量比较关系成立值为1，否则为0；b.两个同维数组比较第28页，讲稿共81页，2023年5月2日，星期三每一对应元素间进行比较，结果为同维的0-1矩阵；c.一标量与数组比较标量与数组中的每一元素作比较，结果为0-1矩阵。A=[35,1,6;26,19,24;3,32,7];B=9;C=9*eye(size(A));C=900090009D=rem(A,3)==0D=001001100D1=C<=AD1=111111110例：D2=B>AD2=011000100第29页，讲稿共81页，2023年5月2日，星期三（3）find函数找出向量或数组中非零元素的位置标识，数组按列顺序找。格式：i=find(x)i为非零元素的序号（以列为顺序）[i,j]=find(x)i,j分别为非零元素的行号和列号A=[3,5,-4,0;5,9,-1,5];i=find(A<=0)[j,k]=find(A<=0)i=567j=121k=334第30页，讲稿共81页，2023年5月2日，星期三2、逻辑运算（2）运算规则&,|,~,&&,||（1）逻辑运算符

a.非零元素的逻辑量为真，其代码为1，否则为0。b.两个标量比较

a&b,真真为真（1），否则为0；a|b,其一真则真（1），否则为0；~a,真变假（0），假变真（1）。&&当运算符左边为1时，才继续执行符号右边的运算||当运算符左边为0时，才继续执行符号右边的运算第31页，讲稿共81页，2023年5月2日，星期三c.两个同维数组比较同一位置按标量规则进行运算，结果为同维的0-1数组。d.标量与数组标量与数组中每个元素作逻辑运算，结果为的0-1数组。e.逻辑运算符的优先级~&

|

&&||优先级依次递减，（3）逻辑函数all(A):A为向量，若其全为非零元素，值为1，否则为0；A为数组，则对A每一列以向量规则运算，值为0-1向量。

第32页，讲稿共81页，2023年5月2日，星期三数学运算、关系运算和逻辑运算三者优先级顺序(P62)：‘,^,.’,.^~,*,/,\,.*,./,.\,+,-,:<,<=,>,>=,~=,==&,|,&&,||any(A):A为向量，若有一元素非零，值为1，否则为0；A为数组，则对A每一列以向量规则运算，值为0-1向量。isinf(A):返回一个与A同维的数组,中相应位置元素的值为无穷大时的值为1,否则为0.isnan(A):返回一个与A同维的数组,

isfinite(A):返回一个与A同维的数组,

第33页，讲稿共81页，2023年5月2日，星期三六、线性代数方程组求解matlab中有两种除运算左除和右除。对于方程Ax=b，A为m×n矩阵，有三种情况：当m=n时，此方程成为“恰定”方程当m>n时，此方程成为“超定”方程当m<n时，此方程成为“欠定”方程

matlab定义的除运算可以很方便地解上述三种方程第34页，讲稿共81页，2023年5月2日，星期三对于方程Ax=0，称为齐次方程，其有两种可能的解：（a）零解，充要条件是R(A)=n,若是方阵,则|A|0;（b）非零解，充要条件是R(A)<n(有基础解系)齐次方程Ax=0的求解利用行阶梯函数:R=rref(A)或有理基函数Z=null(A,‘r’)(即基础解系)r表示Z的列向量是方程Ax=0的有理基.1.齐次方程组的求解第35页，讲稿共81页，2023年5月2日，星期三A=[1,1,2,-1;2,2,4,-2;2,1,1,-1;2,2,1,2];formatratR=rref(A)例:

x1+x2+2x3-x4=02x1+2x2+4x3-2x4=02x1+x2+x3-x4=0

2x1+2x2+x3+2x4=0结果：R=

100-4/30103001-4/30000

第36页，讲稿共81页，2023年5月2日，星期三通解Z=null(A,'r')结果：Z=4/3-34/31第37页，讲稿共81页，2023年5月2日，星期三2.非齐次方程组的求解非齐次线性方程组求解步骤:（1）判断Ax=b是否有解,若有则继续二步;（2）求Ax=b的一个特解;若是恰定方程,x=inv(A)b(det(A)~=0)或x=A\b运算求解方程,而超定方程和欠定方程x=A\b

;（3）求Ax=0的通解;求基础解系:Z=null(A,‘r’)（4）Ax=b的通解等于Ax=0的通解与Ax=b的一个特解相加。对于方程Ax=b，称为非齐次方程。Ax=b其解存在三种可能：（a）无解，充要条件是R(A)<R(A,b)（b）唯一解，充要条件是R(A)=R(A,b)=n（c）无穷解，充要条件是R(A)=R(A,b)<n第38页，讲稿共81页，2023年5月2日，星期三例:

x1+x2-3x3-x4=12x1+2x2+x3-2x4=32x1+x2+x3-x4=2

A=[1,1,-3,-1;2,2,1,-2;2,1,1,-1];b=[1;3;2];B=[A,b];ra=rank(A);rb=rank(B);formatratifra==rb&ra==length(A(1,:))x=A\belseifra==rb&ra<length(A(1,:))x=A\b,c=null(A,'r')else'Theequationhasnosolve.'endx=3/711/70

c=0101结果：通解：第39页，讲稿共81页，2023年5月2日，星期三解线性方程组的一般函数文件如下：function[x,y]=line_solution(A,b)[m,n]=size(A);y=[];ifnorm(b)>0%非齐次方程组

ifrank(A)==rank([A,b])%方程组相容

ifrank(A)==n%有唯一解

x=A\b;

else%方程组有无穷多个解,基础解系

disp('原方程组有有无穷个解，其齐次方程组的基础解系为y，特解为x');y=null(A,'r');x=A\b;

end第40页，讲稿共81页，2023年5月2日，星期三

else%方程组不相容，无解

disp(‘方程组无解');end

else%齐次方程组

ifrank(A)>=n%列满秩

x=zero(m,1)%0解

else%非0解

disp('方程组有无穷个解，基础解系为x');x=null(A,'r');endendreturn第41页，讲稿共81页，2023年5月2日，星期三3.非线性方程的解法二分法：erfen(‘fun’,x0,xf,esp)迭代法：newtonger4.非线性方程组的解法不动点迭代法：iterateproNewton迭代法：newtonpro第42页，讲稿共81页，2023年5月2日，星期三1、多项式创建

f(x)=anxn+an-1xn-1+……+a0

可用行向量p=[anan-1……a1+a0]表示六、多项式运算（P111）

2、多项式运算（1）roots——求多项式的根例：求x3+4x2+2x+1的根p=[1,4,2,1];r=roots(p)第43页，讲稿共81页，2023年5月2日，星期三结果：r=-3.5115-0.2442+0.4745i-0.2442-0.4745i可利用poly指令和r返回多项式形式pp=poly(r);f=poly2sym(pp)或f=poly2str(pp,’x’)—函数文件，显示数学多项式的形式结果：x^3+4*x^2+2*x+1第44页，讲稿共81页，2023年5月2日，星期三（2）求多项式的值

格式：f=polyval(p,s)p为多项式系数向量，s为数值或数组（3）多项式乘、除运算(向量的卷积和解卷积)乘（卷积）：c=conv(a,b)[d,r]=deconv(c,a)余数c除a后的整数除（解卷积）：第45页，讲稿共81页，2023年5月2日，星期三n=1:50;%定义序列的长度hb=zeros(1,50);hb(1)=1;hb(2)=2.5;hb(3)=2.5;hb(4)=1;closeall;subplot(3,1,1);stem(hb);title(‘系统hb[n]’);m=1:50;%定义序列的长度例：第46页，讲稿共81页，2023年5月2日，星期三A=444.128;%设置信号的有关参数a=50*sqrt(2.0)*pi;T=0.001;%采样率w0=50*sqrt(2.0)*pi;x=A*exp(-a*m*T).*sin(w0*m*T);subplot(3,1,2);stem(x);title(‘输入信号x[n’);y=conv(x,hb);subplot(3,1,3);stem(y);title(‘输出信号y[n]’);第47页，讲稿共81页，2023年5月2日，星期三第48页，讲稿共81页，2023年5月2日，星期三（4）多项式微分命令格式：polyder(p):求p的微分polyder(a,b):求多项式a(x)*b(x)乘积的微分[p,q]=polyder(a,b):求多项式a(x)/b(x)商的微分例：a=[12345];poly2str(a,'x')ans=x^4+2x^3+3x^2+4x+5b=polyder(a)结果：b=4664poly2str(b,'x')结果：ans=4x^3+6x^2+6x+4第49页，讲稿共81页，2023年5月2日，星期三七、数值微积分及微分方程求解1、差分与微分dx=diff(x,k)结果为(m-k)xn维k阶差分阵dyx=diff(y)./diff(x)2、积分trapz(x,y)梯形法积分quad(‘函数文件’,x0,xf,tol,trace)

自适应递推辛普生法（低阶）quad8(‘函数文件’,x0,xf,tol,trace)自适应牛顿—柯西法（高阶）dblquad8(‘函数文件’,x0,xf,y0,yf,

tol,trace)二重定积分第50页，讲稿共81页，2023年5月2日，星期三formatshorts=quad8('jifen',0,pi)functiony=jifen(x)y=x.*sin(x)./(1+(cos(x)).^2);s=2.4674求：取jifen文件名保存另开编辑窗口第51页，讲稿共81页，2023年5月2日，星期三微分方程求解的仿真算法有多种，常用的有Euler（欧拉法）、RungeKutta(龙格-库塔法。Euler法称一步法，用于一阶微分方程3、微分方程求解格式：[x,y]=eulerpro(‘fun’,x0,xf,y0,h)第52页，讲稿共81页，2023年5月2日，星期三龙格-库塔法：实际上取两点斜率的平均斜率来计算的，其精度高于欧拉算法。龙格-库塔法：ode23ode45[x,y]=ode23(‘函数文件’,x0,xf,y0,tol,trace)[x,y]=ode45(‘函数文件’,x0,xf,y0,tol,trace)说明：（1）函数文件的书写中，必须书写成一阶微分方程dy/dx=f(x,y)的形式;函数文件须以file(x,y)形式；（2）微积分的上下限，是已知的初始条件，其为列向量，且高阶项在前。第53页，讲稿共81页，2023年5月2日，星期三例：求x’’+(x2-1)x’+x=0,已知x|t=0=0，x’|t=0=0.25.为方便令x1=x，x2=x分别对x1,x2求一阶导数，整理后写成一阶微分方程组形式

x’1=x2x’2=x2(1-x12)-x1建立m文件解微分方程第54页，讲稿共81页，2023年5月2日，星期三建立m文件functionxdot=wf(t,x)xdot=zeros(2,1)xdot(1)=x(2)xdot(2)=x(2).*(1-x(1).^2)-x(1)给定区间、初始值;求解微分方程t0=0;tf=20;x0=[0.25,0]';[t,x]=ode23('wf',t0,tf,x0)plot(t,x)第55页，讲稿共81页，2023年5月2日，星期三第56页，讲稿共81页，2023年5月2日，星期三八、拟合与插值多项式拟合利用最小二乘法，根据已知数据拟合出多项式系数格式：p=polyfit(x,y,n)例：已知五组数据(1,5),(2,43),(3,128),(4,280),(5,500),试拟合数据的趋势x=1:5;y=[5,43,128,280,500];p=polyfit(x,y,3)结果：p=1.750015.0357-20.71439.2000第57页，讲稿共81页，2023年5月2日，星期三x1=0:0.1:6;y1=polyval(p,x1);plot(x,y,'*',x1,y1,'-r')第58页，讲稿共81页，2023年5月2日，星期三2.插值插值的定义——是对某些集合给定的数据点之间函数的估值方法。当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具。Matlab提供了一维、二维、三次样条等许多插值选择第59页，讲稿共81页，2023年5月2日，星期三（1）Lagrange插值

对于(x1,y1)，…，(xn,yn)已知值，可求出x1~xn间任何x对应下的y值。格式：y=lagrange(x0,y0,x)第60页，讲稿共81页，2023年5月2日，星期三利用已知点确定未知点粗糙——精确集合大的——简化的（2）分段线性插值

intep1、interp2、linear、spline三次样条插值、cubic三次插值。（3）Hermite插值

要求节点上函数值相等，导数值相等，甚至要求高阶导数值相等。格式：y=hermite(x0,y0,y1,x)

y1指已知点的一阶导数值。第61页，讲稿共81页，2023年5月2日，星期三八、数据分析与优化问题1、数据分析（1）常用指令max(A)——求向量或矩阵各列最大值

[y,m]=max(A)—y记录各列最大值

,m记录各列最大值

的行号

[y,m]=

max(A,[],dim)—dim为1结果同上

,为2则记录各行最大值与相应的列号min(A)——求向量或矩阵各列最小值第62页，讲稿共81页，2023年5月2日，星期三mean(A)——求向量或矩阵各列的平均值

mean(A,dim)——dim的值和意义同上所示median(A)——求向量或矩阵各列的中值

median(A,dim)——dim的值和意义同上所示sum(A)——求向量或矩阵各列的和

sum(A,dim)[y,m]=min(A)

[y,m]=

min(A,[],dim)——

y,m和dim的值和意义同上所示第63页，讲稿共81页，2023年5月2日，星期三std——各列标准方差

std(A,flag)std(A,flag,dim)prod——各列求积,如prod(1:10)即为10！

prod(A,dim)sort(A)

——

各列递增排序，

sort(A,dim,mode)——

dim的值取1和2分别对应按列排或按行排;mode为排序的方式,若取’ascend’为升序,取为’descend’降序.第64页，讲稿共81页，2023年5月2日，星期三cov——协方差（每列间求方差）常用corrcoef——相关阵（每列间求相关系数）常用第65页，讲稿共81页，2023年5月2日，星期三（2）回归分析目的：根据实验或检测的数据，建立回归方程，找出应变量x1,

x2….,xn与响应值y1,

y2….,yn间的相关关系，来检测数据的可靠性。第66页，讲稿共81页，2023年5月2日，星期三步骤：回归方程系数的确定利用最小二乘法确定，常用polyfit(x,y,n)b.

可靠性检验

相关系数法：R=corrcoef(x,y),r=R(1,2)若则可靠，相关性显著。此法适用于一元线性回归。第67页，讲稿共81页，2023年5月2日，星期三F检验法：(不介绍）例：电容器充电时，电压达到100伏后而放电，测得t与u的关系：t(s)012345678910u(V)100755540302015101055试建立与的回归方程，并判别其可靠性。%由画图可知u=aebt,lnu=lna+btt=0:10;u=[100,75,55,40,30,20,15,10,10,5,5];p=polyfit(t,log(u),1)第68页，讲稿共81页，2023年5月2日，星期三b=p(1),a=exp(p(2)),pr=corrcoef(t,log(u));r=pr(1,2),结果：

r=-0.9995

%显示回归方程disp(['u=',num2str(a),‘*exp(',num2str(b),…‘*t)’])结果：u=100.7890*exp(-0.31264*t)%画回归图tt=0:0,1:13;uu=exp(polyval(p,tt));plot(t,u,’*’,tt,uu)第69页，讲稿共81页，2023年5月2日，星期三第70页，讲稿共81页，2023年5月2日，星期三2、最优化问题