从零开始学MATLAB（六）数值微积分与方程求解

从零开始学MATLAB（六）数值微积分与方程求解

1.数值微分与数值积分

（1）数值微分
MATLAB提供了求向前差分的函数diff，其调用格式有三种：
①dx=diff(x)：计算向量x的一阶向前差分，dx(i)=x(i+1)-x(i),i=1,2,…,n-1
②dx=diff(x,n)：计算向量x的n阶向前差分。
③dx=diff(A,n,dim)：计算矩阵A的n阶差分，dim=1时（默认状态），按列计算差分；dim=2，按行计算差分。
差商便是该点导数的近似值。
（2）数值积分
①基于自适应辛普森方法
[l,n]=quad(filename,a,b,tol,trace)
②基于自适应Gauss-Lobatto方法
[l,n]=quadl(filename,a,b,tol,trace)
其中，filename是被积函数名；a和b分别是定积分的下限和上限，积分限[a,b]必须是有限的，不能为无穷大（inf）；tol用来控制积分精度（默认取10的-6次方）；trace控制是否展现积分过程，取非0则展现，取0则不展现（默认取0）；返回参数l即定积分的值，n为被积函数的调用次数。
③基于全局自适应积分方法：
l=integral(filename,a,b)，其中，l是计算得到的积分；filename是被积函数；a和b分别是定积分的下限和上限，积分限可以为无穷大。
④基于自适应高斯-克朗罗德方法：
[l,err]=quadgk(filename,a,b)，其中，err返回近似误差范围，其他参数的含义和用法与quad函数相同。积分上下限可以是无穷大（-inf或inf），也可以是复数。如果积分上下限是复数，则quadgk函数在复平面上求积分。
⑤基于梯形积分法：l=trapz(x,y)，其中，向量x、y定义函数关系y=f(x)

2.线性方程组求解：直接法、迭代法

（1）直接法：高斯消去法、列主元消去法、矩阵的三角分解法
①利用左除运算符的直接解法
Ax=b---------x=A\b
注：如果矩阵A是奇异的或接近奇异的，则MATLAB会给出警告信息。
②利用矩阵分解求解线性方程组：LU分解、QR分解、Cholesky 分解
LU分解思想：将矩阵表示成一个下三角矩阵和一个上三角矩阵的乘积

LU分解函数：
[L,U]=lu(A)：产生一个上三角振U和一个变换形式的下三角阵L，使之蛮子A=LU。注意，这里的矩阵A必须是方阵。
[L,U,P]=lu(A)：产生一个上三角阵U和一个下三角阵L以及一个置换矩阵P，使之满足PA=LU。同样，矩阵A必须为方阵。

（2）迭代解法：雅克比迭代法、高斯-赛德尔迭代法

雅可比迭代法的函数文件jacobi.m：
function [y,n]=jacobi(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=D(L+U);
f=D\b;
y=Bx0+f;
n=1;
while norm(y-x0)>=ep
x0=y;
y=Bx0+f;
n=n+1;
end

Gauss-Serdel迭代法的函数文件gauseidel.m：
function [y,n]=gauseidel(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1);
B=(D-L)\U;
f=(D-L)\b;
y=Bx0+f;
n=1;
while norm(y-x0)>=ep
x0=y;
y=Bx0+f;
n=n+1;
end

3.非线性方程组求解

（1）单变量非线性方程求解：x=fzero(filename,x0)
其中，filename是待求根方程左端的函数表达式，x0是初始值。
（2）非线性方程组的求解
函数的调用格式：x=fsolve(filename,x0,option)
其中，x为返回的近似解，filename是待求根方程的左端的函数表达式，x0是初值，option用于设置优化工具箱的优化参数，可以调用optimset函数来完成。
（3）函数极值的计算
只提供极小值的求法，极大值通过-f(x)的极小值求。
①无约束最优化问题
[xmin,fmin]=fminbnd(filename,x1,x2,option)
[xmin,fmin]=fminsearch(filename,x0,option)
[xmin,fmin]=fminunc(filename,x0,option)
其中，filename是定义的目标函数。第一个函数的输入变量x1、x2分别表示被研究区间的左、右边界。后两个函数的输入变量x0是一个向量，表示极值点的初值。Option为优化参数，可以通过optimset函数来设置。
②有约束最优化问题
约束条件可细化为：线性不等式约束、线性等式约束、非线性不等式约束、非线性等式约束、x的上界与下界
求有约束条件下最小值的函数为：[xmin,fmin]=fmincon(filename,x0,A,b,Aeq,beq,Lbnd,Ubnd,NonF,option)
其中，xmin、fmin、filename、x0和option的含义与求最小值函数相同。其余参数为约束条件，包括线性不等式约束、线性等式约束、x的下界和上界以及定义非线性约束的函数。如果某个约束不存在，则用空矩阵来表示。

4.常微分方程数值求解

（1）一般概念：单步法（代表是龙格-库塔法）、多步法（代表是亚当斯法）
（2）求解函数：
①[t,y]=solver(filename,tspan,y0,option)
其中，t和y分别给出时间向量和相应的数值解；solver为求常微分方程数值解的函数；filename是定义f(t,y)的函数名，该函数必须返回一个列向量；tspan形式为[t0,tf]，表示求解区间；y0是初始状态向量；option是可选参数，用于设置求解属性。
②常微分方程数值求解函数的统一命名格式：odennxx
其中，Ode是常微分方程，nn是数字，代表所用方法的阶数（前一个数字是描述公式的阶数，后一个字母描述误差估计调节步长的公式阶数）；xx是字母，用于标注方法的专门特征（s表示stiff，刚性问题）。