数值计算功能

1、关于数值计算功能第一张，PPT共八十一页，创作于2022年6月数值运算的功能创建矩阵矩阵运算多项式运算线性方程组数值统计线性插值函数优化微分方程的数值解第二张，PPT共八十一页，创作于2022年6月一、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第三张，PPT共八十一页，创作于20

2、22年6月以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)等第四张，PPT共八十一页，创作于2022年6月二、命令行的基本操作创建矩阵的

3、方法(1) 直接输入法规则：矩阵元素必须用 括住矩阵元素必须用逗号或空格分隔 在 内矩阵的行与行之间必须 用分号分隔第五张，PPT共八十一页，创作于2022年6月 矩阵元素可以是任何matlab表达式 ，可以是实数 ，也可以是复数，复数可用特殊函数i，j 输入 a=1, 2 ,3;4, 5, 6 矩阵元素a = 1 2 3 4 5 6x=2 pi/2;sqrt(3) 3+5i第六张，PPT共八十一页，创作于2022年6月符号的作用逗号和分号的作用 逗号和分号可作为指令间的分隔符，matlab允许多条语句在同一行出现。 分号如果出现在指令后，屏幕上将不显示结果，但存储在工作空间中。当一个指令或矩

4、阵太长时，可用续行第七张，PPT共八十一页，创作于2022年6月 用于选出矩阵指定行、列及元素。a=A(:,2:4)或 a=A(:,2,3,4)冒号的作用：用于生成等间隔的向量，默认间隔为1。i=3:2:7循环语句for n=1:10第八张，PPT共八十一页，创作于2022年6月(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) )

5、 生成ab之间n个数值线性分布的向量logspace(a,b,n) ) 生成10a10b之间按对数等分的n个元素的向量第九张，PPT共八十一页，创作于2022年6月r=eye(3) r=1 0 0 0 1 0 0 0 1r=eye(3,4)r = 1 0 0 0 0 1 0 0 0 0 1 0r=eye(4,3)r = 1 0 0 0 1 0 0 0 1 0 0 0第十张，PPT共八十一页，创作于2022年6月 还有伴随矩阵、稀疏矩阵、魔方矩阵、对角矩阵、范德蒙等矩阵的创建，就不一一介绍了。注意：matlab严格区分大小写字母，因此a与A是两个不同的变量。 第十一张，PPT共八十一页，创作于2

6、022年6月(3) 用M文件创建矩阵适用于较大的矩阵且需经常调用的矩阵。% mymatrix.m creation of matrix JZJZ=1,2,3,4,5;6,7,8,9,10;11,12,13,14,15;取mymatrix名保存若JZ已由其它文件运行后生成而需经常调用，可用save指令保存。save mymatrix JZ(4) 从外部数据文件调入矩阵load mymatrix JZ第十二张，PPT共八十一页，创作于2022年6月三、 矩阵的修改与操作 直接修改 可用键找到所要修改的矩阵，用键移动到要修改的矩阵元素上即可修改。 指令修改 可以用A(,)= 来修改。(1) 矩阵的修

7、改(5) 多维矩阵第十三张，PPT共八十一页，创作于2022年6月例如a=1 2 0;3 0 5;7 8 9a =1 2 0 3 0 5 7 8 9a(3,3)=0a =1 2 0 3 0 5 7 8 0还可以用函数subs、find函数修改。第十四张，PPT共八十一页，创作于2022年6月(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)等第十五张，PPT共八十一页，创作于2022年6

8、月(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) 结果相同 2 3 6 7 10 11第十六张，PPT共八十一页，创作于2022年6月矩阵的扩展a. 利用矩阵标识块赋值指令扩展X(m1:m2,n1:n2)=A生成新阵X，矩阵除了赋值阵A和已存在的元素外，其余为0 b. 利用方括号和小矩阵生成大矩阵例：利用A=1,2,7,8;-1,19,7,10生成44的B阵，使其

9、第1、2行，第1至3列元素取A阵第2列至第4列元素，第3、4行元素全为1，其余为0。第十七张，PPT共八十一页，创作于2022年6月法一：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) 2 7 8 0 19 7 10 0 1 1 1 1 1 1 1 1或B(1:2,1:3)=A(:,2:4);B(3:4,1:4)=ones(2,4)第十八张，PPT共八十一页，创作于2022年6月矩阵的结构变换P91-94rot90(A,k):逆时针旋转fliplr:左右翻flipud:上

10、下翻diag:抽取主对角线或生成对角阵tril(A,k): 抽取主下三角(k=0,+1,-1,)triu(A,k):抽取主上三角(k=0,+1,-1,)B=reshape(A, m, n）根据A阵重组为mn的B阵第十九张，PPT共八十一页，创作于2022年6月矩阵加、减（,）运算规则：与线性代数运算规则相同四、矩阵运算2. 矩阵乘（）运算规则：与线性代数运算规则相同第二十张，PPT共八十一页，创作于2022年6月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.0000 5.0000 6.1416 4.0000 + 1.0000i

11、11.0000 14.0000 10.0000 + 1.0000i 13.0000 a=1 2 3;4 5 6;7 8 0;b=1;2;3;c=a*bc =14 32 23第二十一张，PPT共八十一页，创作于2022年6月 a p a 自乘p次幂方阵1的整数3. 矩阵乘方 an,ap,pa对于p的其它值,计算将涉及特征值和特征向量，如果p是矩阵，a是标量ap使用特征值和特征向量自乘到p次幂；如a,p都是矩阵，ap则无意义。第二十二张，PPT共八十一页，创作于2022年6月inv 矩阵求逆rank 求矩阵的秩det 行列式的值eig 矩阵的特征值diag 对角矩阵 矩阵转置sqrt 矩阵开方 4

12、. 矩阵的其它运算 (P96)第二十三张，PPT共八十一页，创作于2022年6月 数组运算指元素对元素的算术运算，与通常意义上的由符号表示的线性代数矩阵运算不同 数组加减(.+, .-) a.+b a.- b5. 矩阵的数组运算 对应元素相加减（与矩阵加减等效）第二十四张，PPT共八十一页，创作于2022年6月2. 数组乘除(，./，.）ab a，b两数组必须有相同的行 和列两数组相应元素相乘。a=1 2 3;4 5 6;7 8 9;b=2 4 6;1 3 5;7 9 10;a.*bans = 2 8 18 4 15 30 49 72 90 第二十五张，PPT共八十一页，创作于2022年6月a

13、./b=b.a a的元素被b的对应元素除a.b=b./a b的元素被a的对应元素除例: a=1 2 3;b=4 5 6; c1=a.b; c2=b./ac1 = 4.0000 2.5000 2.0000c2 = 4.0000 2.5000 2.0000第二十六张，PPT共八十一页，创作于2022年6月3. 数组乘方(.) 元素对元素的幂例:a=1 2 3;b=4 5 6;z=a.2z = 1.00 4.00 9.00z=a.bz = 1.00 32.00 729.00第二十七张，PPT共八十一页，创作于2022年6月五、关系运算、逻辑运算及其函数 关系运算和逻辑运算均按照数组运算的规则和定义进

14、行的。1、关系运算（2）运算规则, =, =, =, =（1）关系运算符 a.两个标量比较 关系成立值为1，否则为0； b.两个同维数组比较第二十八张，PPT共八十一页，创作于2022年6月每一对应元素间进行比较，结果为同维的0-1矩阵；c. 一标量与数组比较标量与数组中的每一元素作比较，结果为0-1矩阵。A=35,1,6;26,19,24;3,32,7;B=9;C=9*eye(size(A);C=9 0 0 0 9 0 0 0 9D=rem(A,3)=0D= 0 0 1 0 0 1 1 0 0D1=CAD2=0 1 1 0 0 0 1 0 0第二十九张，PPT共八十一页，创作于2022年6月

15、（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 = 5 6 7j = 1 2 1k = 3 3 4第三十张，PPT共八十一页，创作于2022年6月2、逻辑运算（2）运算规则&, |, , &, |（1）逻辑运算符 a.非零元素的逻辑量为真，其代码为1，否则为0。 b.两个标量比较 a&b, 真真为真（1），否则为0； a|b, 其一真则真（1），否则为0； a,

16、真变假（0），假变真（1）。&当运算符左边为1时，才继续执行符号右边的运算|当运算符左边为0时，才继续执行符号右边的运算第三十一张，PPT共八十一页，创作于2022年6月c.两个同维数组比较同一位置按标量规则进行运算，结果为同维的0-1数组。 d.标量与数组 标量与数组中每个元素作逻辑运算，结果为的0-1数组。 e.逻辑运算符的优先级& | & |优先级依次递减，（3）逻辑函数all(A):A为向量，若其全为非零元素，值为1，否则 为0；A为数组，则对A每一列以向量规则运算，值为0-1向量。 第三十二张，PPT共八十一页，创作于2022年6月数学运算、关系运算和逻辑运算三者优先级顺序(P62)

17、：, , ., . , *, /, , .*, ./, ., +, -, : , , =, =, = &, |, &, |any(A):A为向量，若有一元素非零，值为1，否则 为0；A为数组，则对A每一列以向量规则运算，值为0-1向量。isinf(A):返回一个与A同维的数组, 中相应位置元素的值为无穷大时的值为1, 否则为0.isnan(A):返回一个与A同维的数组, isfinite(A):返回一个与A同维的数组, 第三十三张，PPT共八十一页，创作于2022年6月六、线性代数方程组求解matlab中有两种除运算左除和右除。对于方程Ax=b，A 为mn矩阵，有三种情况： 当m=n时，此方程

18、成为“恰定”方程 当mn时，此方程成为“超定”方程 当mn时，此方程成为“欠定”方程 matlab定义的除运算可以很方便地解上述三种方程第三十四张，PPT共八十一页，创作于2022年6月对于方程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. 齐次方程组的求解第三十五张，PPT共八十一页，创作于2022年6月A=1,1,2,-1;2

19、,2,4,-2;2,1,1,-1;2,2,1,2;format ratR=rref(A)例: x1+x2+2x3-x4=0 2x1+2x2+4x3-2x4=0 2x1+x2+x3-x4=0 2x1+2x2+x3+2x4=0结果：R= 1 0 0 -4/3 0 1 0 3 0 0 1 -4/3 0 0 0 0 第三十六张，PPT共八十一页，创作于2022年6月通解Z=null(A,r)结果：Z= 4/3 -3 4/3 1第三十七张，PPT共八十一页，创作于2022年6月2. 非齐次方程组的求解非齐次线性方程组求解步骤:（1）判断Ax=b是否有解, 若有则继续二步;（2） 求Ax=b的一个特解;若

20、是恰定方程, x=inv(A)b (det(A)=0)或x=Ab 运算求解方程, 而超定方程和欠定方程x=Ab ; （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第三十八张，PPT共八十一页，创作于2022年6月例: x1+x2-3x3-x4=1 2x1+2x2+x3-2x4=3 2x1+x2+x3-x4=2

21、 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); format ratif ra=rb&ra=length(A(1,:) x=Abelseif ra=rb&ra0 %非齐次方程组 if rank(A)=rank(A,b) %方程组相容 if rank(A)=n %有唯一解 x=Ab; else %方程组有无穷多个解,基础解系 disp(原方程组有有无穷个解，其齐次方程组的基础解系为y，特解为x); y=null(A,r); x=Ab; end第四十张，PPT共八十一页，创作于2022年6月 else %方程组

22、不相容，无解 disp(方程组无解); end else %齐次方程组 if rank(A)=n %列满秩 x=zero(m,1) %0解 else %非0解 disp(方程组有无穷个解，基础解系为x); x=null(A,r); end endreturn第四十一张，PPT共八十一页，创作于2022年6月3.非线性方程的解法二分法：erfen(fun,x0,xf,esp)迭代法：newton ger4.非线性方程组的解法不动点迭代法：iterateproNewton迭代法：newtonpro第四十二张，PPT共八十一页，创作于2022年6月1、多项式创建 f(x)=anxn+an-1xn-1

23、+a0 可用行向量 p=an an-1 a1 +a0表示六、 多项式运算（P111） 2、多项式运算（1）roots 求多项式的根例：求x3+4x2+2x+1的根p=1,4,2,1;r=roots(p)第四十三张，PPT共八十一页，创作于2022年6月结果：r = -3.5115 -0.2442 + 0.4745i -0.2442 - 0.4745i可利用poly指令和r返回多项式形式pp=poly(r); f=poly2sym(pp)或f=poly2str(pp,x) 函数文件，显示数学多项式的形式结果：x3+4*x2+2*x+1第四十四张，PPT共八十一页，创作于2022年6月（2）求多项

24、式的值 格式：f=polyval(p,s)p为多项式系数向量，s为数值或数组（3） 多项式乘、除运算(向量的卷积和解卷积)乘（卷积）：c=conv (a,b)d,r=deconv(c,a)余数c除a后的整数除（解卷积）：第四十五张，PPT共八十一页，创作于2022年6月 n=1:50; %定义序列的长度 hb=zeros(1,50); hb(1)=1;hb(2)=2.5;hb(3)=2.5;hb(4)=1; close all; subplot(3,1,1);stem(hb);title(系统 hbn );m=1:50; %定义序列的长度例：第四十六张，PPT共八十一页，创作于2022年6月A

25、=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(输入信号 xn);y=conv(x,hb); subplot(3,1,3);stem(y);title(输出信号yn);第四十七张，PPT共八十一页，创作于2022年6月第四十八张，PPT共八十一页，创作于2022年6月（4）多项式微分命令格式：polyder(p): 求p的微分polyder(a,b): 求多项式a(x)*b(x)

26、乘积的微分p,q=polyder(a,b): 求多项式a(x)/b(x)商的微分例：a=1 2 3 4 5; poly2str(a,x)ans = x4 + 2 x3 + 3 x2 + 4 x + 5b=polyder(a) 结果：b = 4 6 6 4poly2str(b,x)结果：ans =4 x3 + 6 x2 + 6 x + 4第四十九张，PPT共八十一页，创作于2022年6月七、数值微积分及微分方程求解1、差分与微分dx=diff(x,k) 结果为(m-k)xn维k阶差分阵dyx=diff(y). /diff(x)2、积分trapz(x,y) 梯形法积分quad( 函数文件,x0,x

27、f, tol,trace) 自适应递推辛普生法（低阶）quad8 ( 函数文件,x0,xf, tol,trace) 自适应牛顿柯西法（高阶）dblquad8 ( 函数文件,x0,xf, y0,yf, tol,trace) 二重定积分第五十张，PPT共八十一页，创作于2022年6月format shorts=quad8(jifen,0,pi)function y=jifen(x)y=x.*sin(x)./(1+(cos(x).2);s = 2.4674求：取jifen文件名保存另开编辑窗口第五十一张，PPT共八十一页，创作于2022年6月微分方程求解的仿真算法有多种，常用的有Euler（欧拉法）

28、、Runge Kutta(龙 格-库塔法。Euler法称一步法，用于一阶微分方程3、微分方程求解格式：x,y=eulerpro(fun,x0,xf,y0,h)第五十二张，PPT共八十一页，创作于2022年6月龙格-库塔法：实际上取两点斜率的平均 斜率来计算的，其精度高于欧拉算法 。龙格-库塔法：ode23 ode45x,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）微积分的上下限，

29、是已知的初始条件，其为列向量，且高阶项在前。第五十三张，PPT共八十一页，创作于2022年6月例：求x+(x2-1)x+x=0, 已知x|t=0=0，x|t=0=0.25.为方便令x1=x，x2=x分别对x1,x2求一阶导数，整理后写成一阶微分方程组形式 x1=x2 x2=x2(1-x12)-x1建立m文件解微分方程第五十四张，PPT共八十一页，创作于2022年6月建立m文件function xdot=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.2

30、5 ,0;t,x=ode23(wf, t0, tf, x0)plot(t,x)第五十五张，PPT共八十一页，创作于2022年6月第五十六张，PPT共八十一页，创作于2022年6月八、拟合与插值多项式拟合利用最小二乘法，根据已知数据拟合出多项式系数格式：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.7500 15.0357 -20.7143 9.2000第五十七张，PPT共八十一页，创作于2022年6月x1

31、=0:0.1:6;y1=polyval(p,x1);plot(x,y,*,x1,y1,-r)第五十八张，PPT共八十一页，创作于2022年6月2.插值插值的定义是对某些集合给定的数据点之间函数的估值方法。当不能很快地求出所需中间点的函数时，插值是一个非常有价值的工具。Matlab提供了一维、二维、 三次样条等许多插值选择第五十九张，PPT共八十一页，创作于2022年6月（1）Lagrange插值 对于(x1,y1)， ，(xn,yn)已知值，可求出x1 xn间任何x对应下的y值。格式：y=lagrange(x0,y0,x)第六十张，PPT共八十一页，创作于2022年6月利用已知点确定未知点粗糙

32、 精确集合大的 简化的（2）分段线性插值 intep1 、interp2、linear、spline三次样条插值、cubic三次插值。（3）Hermite插值 要求节点上函数值相等，导数值相等，甚至要求高阶导数值相等。格式：y=hermite(x0,y0,y1,x) y1指已知点的一阶导数值。第六十一张，PPT共八十一页，创作于2022年6月八、数据分析与优化问题1、数据分析（1）常用指令max(A) 求向量或矩阵各列最大值 y, m=max(A) y记录各列最大值 , m记录 各列最大值 的行号 y, m= max(A, , dim) dim为1结果同上 , 为2则记录各行最大值 与相应的列

33、号min (A) 求向量或矩阵各列最小值第六十二张，PPT共八十一页，创作于2022年6月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的值和意义同上所示第六十三张，PPT共八十一页，创作于2022年6月std 各列标准方差 std(A, flag) std(A, flag, dim)prod 各列求积

34、,如prod(1:10)即为10！ prod (A , dim)sort(A) 各列递增排序， sort(A, dim, mode) dim的值取1和2分别对应按列排或按行排; mode为排序的方式, 若取ascend为升序, 取为descend降序.第六十四张，PPT共八十一页，创作于2022年6月cov 协方差（每列间求方差）常用corrcoef 相关阵（每列间求相关系数）常用第六十五张，PPT共八十一页，创作于2022年6月（2）回归分析目的：根据实验或检测的数据，建立回归方程，找出应变量x1 , x2 .,xn与响应值y1 , y2 .,yn间的相关关系，来检测数据的可靠性。第六十六张

35、，PPT共八十一页，创作于2022年6月步骤：回归方程系数的确定利用最小二乘法确定，常用polyfit(x,y,n)b. 可靠性检验 相关系数法 ：R=corrcoef(x,y), r=R(1,2)若则可靠，相关性显著。此法适用于一元线性回归。第六十七张，PPT共八十一页，创作于2022年6月F检验法 ：(不介绍）例：电容器充电时，电压达到100伏后而放电，测得t与u的关系：t(s)012345678910u(V)100755540302015101055试建立与的回归方程，并判别其可靠性。%由画图可知u=aebt, lnu=lna+btt=0:10;u=100,75,55,40,30,20,

36、15,10,10,5,5;p=polyfit(t,log(u),1)第六十八张，PPT共八十一页，创作于2022年6月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)第六十九张，PPT共八十一页，创作于2022年6月第七十张，PPT共八十一页，创作于2022年6月2、最优化问题 在工程设计、经济管理和科学研究等许多领域中，需要在一切可能的方案中选择最好最优方案；在安排生产或销售中，如何安排人力或物力，才能使产值最高或利润最高、库存成本最小；在投资计划中，如何投资收益最高等等，均为最优化问题。 最优化问题一般包括线性优化和非线性优化两类。第七十一张，PPT共八十一页，创作于