一、数值微分与数值积分
1、数值微分
MATLAB提供了求向前差分的函数diff，调用格式：
(1)dx=diff(x)：计算向量x的向前差分，dx(i)=x(i+1)-x(i),i=1,2,……,n-1。
(2)dx=diff(x,n)：计算向量x的n阶向前差分，例如，diff(x,2)=diff(diff(x))。
(3)dx=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时(默认状态)，按列计算差分；dim=2时，按行计算差分。
注意：diff函数计算的是向量元素间的差分，故差分向量元素的个数比原向量少了一个。同样对于矩阵，差分后的矩阵比原矩阵少了一行或者一列。
另外，计算差分之后，可以用f(x)在某点的差商作为其导数的近似值。
例子：设f(x)=sinx，在[0,2Π]范围内随机采样，计算f’(x)的近似值，并与理论f’(x)=cosx进行比较
x=linspace(0,2*pi,5000); y=sin(x); f1=diff(y)./diff(x); f2=cos(x(1:end-1));%差分向量元素比原向量元素少了一个 plot(x(1:end-1),f1,'r',x(1:end-1),f2,'b'); d=norm(f1-f2)%求f1,f2的2范数
2、数值积分
(1)数学原理
常用的求定积分的方式是利用牛顿-莱布尼兹公式：
但是当被积函数的原函数无法用初等函数表示，或者被积函数是用离散的形式给出，这个时候需要用数值解法求解定积分。
求定积分的数值解法很多种，如梯形法、辛普森法、高斯求积法等等。基本思想都是将积分区间分成n个子区间，这样定积分问题就分解成每个子区间上积分后求和的问题：
(2)数值积分的实现
- 基于自适应辛普森方法
[I,n]=quad(filename,a,b,tol,trace)
- 基于自适应Gauss-Lobattofangf
[I,n]=quadl(filename,a,b,tol,trace)
其中，filename是被积函数名；a和b分别是定积分的下限和上限，积分限[a,b]必须是有限的，不能为无穷大；tol用来控制积分精度，默认取
；trace控制是否展现积分过程，若取非0则展现，取0则不展现，默认取0；返回参数I即定积分的值，n为被积函数的调用次数。
例子：分别用quad函数和quadl函数求定积分的近似值，并在相同的积分精度下比较被积函数的调用次数。
format long f=@(x)4./(1+x.^2); [I1,n1]=quad(f,0,1,1e-8) [I2,n2]=quadl(f,0,1,1e-8) (atan(1)-atan(0))*4%输出理论值
基于全局自适应积分方法
I=integral(filename,a,b)
其中，I是计算得到的积分；filename是被积函数；a和b分别是定积分的下限和上限，积分限可以为无穷大。
- 基于自适应的Gauss-Lobattofangf(高斯-克朗罗德)方法
[I,err]=quadgk(filename,a,b)
其中，err为返回近似误差范围，其他参数与quad函数相同。积分上下限可以是无穷大(-Inf或Inf)，也可以是复数。如果积分上下限是复数，则quadgk函数在复平面上求积分。
f=@(x)sin(1./x)./x.^2; I=quadgk(f,2/pi,+Inf)%求得I=1.0000
- 基于梯形积分法
其中，向量x、y定义函数关系y=f(x)
trapz函数采用梯形积分法则，积分的近似值为：
例子：设x=1:6，y=[6,8,11,7,5,2]，用trapz函数计算定积分
x=1:6; y=[6,8,11,7,5,2]; I1=trapz(x,y)%求得I1=35 I2=sum(diff(x).*(y(1:end-1)+y(2:end))/2)%求得I2=35
(3)多重定积分的数值求解
- 求二重积分的数值解：
I=integral2(filename,a,b,c,d)
I=quad2d(filename,a,b,c,d)
I=dbquad(filename,a,b,c,d,tol)
- 求三重积分的数值解：
I=integral3(filename,a,b,c,d,e,f)
I=triplrquad(filename,a,b,c,d,e,f,tol)
f1=@(x,y)exp(-x.^2/2).*sin(x.^2+y); I1=quad2d(f1,-2,2,-1,1)%求得I1=1.5745 f2=@(x,y,z)4*x.*z.*exp(-z.^2.*y-x.^2); I2=integral3(f2,0,pi,0,pi,0,1)%求得I2=1.7328
二、线性方程求解
1、直接法
- 高斯（Gauss）消去法
- 列主元消去法
- 矩阵的三角分解法
(1)利用左除符号
MATLAB提供了一个左除运算符“\”用于求解线性方程组。如线性方程组Ax=b → x=A-1b→x=A\b
如果矩阵A是奇异或者接近奇异的，会给出警告信息。
例子：用左除运算符求解下列线性方程组
(2)利用矩阵分解求解线性方程组
矩阵分解是将一个给定的矩阵分解成若干个特殊类型矩阵的乘积，从而将一个一般的矩阵计算问题转化为几个易求得特殊矩阵的计算问题。此方法优点是运算速度快，可以节省存储空间。
- LU分解
①矩阵的LU分解是将一个n阶矩阵A表示为一个下三角矩阵L和一个上三角矩阵U的乘积。只要方阵是非奇异的，LU分解总可以进行。
Ax=b→Ly=b，Ux=y
②MATLAB提供函数lu，有两种调用格式：
[L,U]=lu(A)：产生一个上三角阵U和一个变换形式的下三角阵L，使之满足A=LU。注意，这里的矩阵A必须是方阵。
[L,U,P]=lu(A)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PA=LU。注意，这里的矩阵A必须是方阵。
③所以，可以利用LU分解求解线性方程组：
Ax=b→LUx=b→x=U(L\b)
或
Ax=b→PAx=Pb→LUx=Pb→x=U(L\P*b)
例子：用LU分解求解下列线性方程组
- QR分解
- Cholesky分解
2、迭代解法
(1)雅可比（Jacobi）迭代法
对于Ax=b，将A分解成（D-L-U），其中D为对角矩阵，L为负的下三角阵，U为负的上三角阵→(D-L-U)x=b
- 求解公式为：
与之对应的迭代公式为：
- 雅可比迭代法的函数文件jacobi.m
function [y,n]=jacobi(A,b,x0,ep)%xo为迭代初值,ep为迭代精度 D=diag(diag(A));%生成A的对角阵； L=-tril(A,-1);%生成A的负的下三角阵 U=-triu(A,1);%生成A的负的上三角阵 B=D\(L+U); f=D\b; y=B*x0+f;%根据初值x0,求第一次迭代 n=1; while norm(y-x0)>=ep%根据前后两次迭代结果的2范数是否很接近来判断 x0=y; y=B*x0+f; n=n+1; end
(2)高斯-赛德尔迭代法(Gauss-Serdel)
- 迭代公式为：
- 雅可比迭代法的函数文件GaussSerdel.m
function [y,n]=jacobi(A,b,x0,ep)%xo为迭代初值,ep为迭代精度 D=diag(diag(A));%生成A的对角阵； L=-tril(A,-1);%生成A的负的下三角阵 U=-triu(A,1);%生成A的负的上三角阵 B=(D-L)\U; f=(D-L)\b; y=B*x0+f;%根据初值x0,求第一次迭代 n=1; while norm(y-x0)>=ep%根据前后两次迭代结果的2范数是否很接近来判断 x0=y; y=B*x0+f; n=n+1; end
(3)例子
- 分别用雅可比迭代法和高斯赛德尔迭代法求解方程组。迭代初值为0，精度为
A=[4,-2,-1;-2,4,3;-1,-3,3]; b=[1,5,0]'; [x,n]=jacobi(A,b,[0,0,0]',1.0e-6) [x,n]=GaussSerdel(A,b,[0,0,0]',1.0e-6)
- 分别用雅可比迭代法和高斯赛德尔迭代法求解方程组。迭代初值为0，精度为
A=[1,2,-2;1,1,1;2,2,1]; b=[9;7;6]; [x,n]=jacobi(A,b,[0;0;0],1.0e-6 [x,n]=GaussSerdel(A,b,[0;0;0],1.0e-6)
- 对比上述两个例子，雅可比迭代法和高斯赛德尔迭代法，其是否收敛、收敛速度的快慢，与实际线性方程组的结构有关
3、直接法与迭代法的对比
- 直接法：以矩阵初等变换为基础，可以求得方程组的精确解；但是占用的空间内存大、程序实现较为复杂；一般适合求解低阶稠密线性方程组。
- 迭代法：从给定初始值逐步逼近精确值的过程，求解过程占用存储空间小，程序设计简单；适合求解大型稀疏矩阵线性方程组；要考虑算法的收敛性。
三、非线性方程求解与函数极值计算
1、非线性方程数值求解
（1）单变量非线性方程求解
x=fzero(filename,x0)
其中，filename是待求根方程左端的函数表达式，x0是初始值。
①例子：
f=@(x)x-1./x+5; x1=fzero(f,-5) x2=fzero(f,1) x3=fzero(f,0.1) fplot(f,[-10,2]) grid on hold on fplot(0,[-10,2]) hold on plot(x1,0,'*',x2,0,'o',x3,0,'p')
所以，利用fzero()函数求解方程，初值的选取很重要。
f=@(x)x.^2-1; x=[]; x0=-0.25:0.001:0.25; for x00=x0 x=[x,fzero(f,x00)]; end plot(x0,x,'-o'); xlabel('初值'); ylabel('方程的根'); axis([-0.25,0.25,-1,1])
结果表明，相同的根对应的处置范围并不连续，求得的根也并非是离初值比较近的根。所以fzero函数执行的是一个数值搜索的过程，搜索结果依赖于函数特性和指定的函数初值。
（2）非线性方程组的求解
x=fsolve(filename,x0,option)
其中，x为返回的近似解，filename是待求根方程左端的函数表达式，x0是初始值，，option用于设置优化工具箱的优化参数，可以调用optimset函数来完成。例如，Display参数设置为’off’时不显示中间结果。
当初值是0.1时，fzero函数无法得到正确的结果，而利用fsolve函数可以。因为不同函数的实现算法不同，适用的场合也不同。
②求下列方程组在(1,1,1)附近的解并对结果进行验证
（因为迭代法是逼近问题，为了展示效果，设置为format long）
2、函数极值的计算
函数极值包括极大值和极小值，或者叫最大值和最小值。MATLAB只考虑吧最小值的计算。
（1）无约束最优化问题
求最小值的函数为：
[xmin,fmin]=fminbnd(filename,x1,x2,option)
[xmin,fmin]=fminsearch(filename,x0,option)
[xmin,fmin]=fminunc(filename,x0,option)
其中，xmin表示极小值点，fmin表示最小值。第一个函数的输入变量x1、x2分别表示被研究区间的左右边界。后两个函数的输入变量x0是一个向量，表示极值点的初值。
（2）有约束最优化问题
即求取一组x，使得目标函数f(x)为最小，且满足约束条件G(x)≤0。
函数为：
[xmin,fmin]=fmincon(filename,x0,A,b,Aeq,beq,Lb,Ub,nonlcon,option)
其中，A,b,Aeq,beq,Lb,Ub,nonlcon分别表示线性不等式约束、线性等式约束、x的下界和上界、定义非线性约束的函数。如果某个约束不存在，用空矩阵表示。
具体各种约束，请查看：https://ww2.mathworks.cn/help/optim/ug/fmincon.html
例子：①求解有约束的最优化问题
f=@(x)0.4*x(2)+x(1)^2+x(2)^2-x(1)*x(2)+1/30*x(1)^3; x0=[0.5;0.5]; A=[-1,-0.5;-0.5,-1];%Ax≤b b=[-0.4;-0.5]; lb=[0;0];%lb≤x≤ub option=optimset('Display','off'); [xmin,fmin]=fmincon(f,x0,A,b,[],[],lb,[],[],option)
例子：②此例特殊说明非线性约束
（x = fmincon(fun,x0,A,b,Aeq,beq,lb,ub,nonlcon) 执行最小化时，满足 nonlcon 所定义的非线性不等式 c(x) 或等式 ceq(x)。fmincon 进行优化，以满足 c(x) ≤ 0 和 ceq(x) = 0。）
编写函数non.m定义非线性约束条件:
function [c,ceq]=non(x); c=[-x(1).^2+x(2)-x(3).^2 x(1)+x(2).^2+x(3).^3-20]; ceq=[-x(1)-x(2).^2+2 x(2)+2*x(3).^2-3];
编写主程序函数
f=@(x)x(1).^2+x(2).^2+x(3).^2+8;; [xmin,fmin]=fmincon(f,rand(3,1),[],[],[],[],zeros(3,1),[],'non')
（3）最小值问题实例
①假设仓库所选点坐标为（x，y），则总里程表达式为：
a=[10,30,16.667,0.555,22.2221]; b=[10,50,29,29.888,49.988]; c=[10,18,20,14,25]; f=@(x)sum(c.*sqrt((x(1)-a).^2+(x(2)-b).^2)); [xmin,fmin]=fminsearch(f,[15,30])
②若地域限制，仓库必须建立在 
上
编写函数non.m定义非线性约束条件:
function [c,ceq]=non(x); c=[]; ceq=[x(2)-x(1)^2];
编写主程序函数
a=[10,30,16.667,0.555,22.2221]; b=[10,50,29,29.888,49.988]; c=[10,18,20,14,25]; f=@(x)sum(c.*sqrt((x(1)-a).^2+(x(2)-b).^2)); option=optimset('Display','off'); [xmin,fmin]=fmincon(f,[15,30],[],[],[],[],[],[],'non',option)
