Matlab数值微积分与方程求解

课程参考：https://www.icourse163.org/learn/CSU-1002475002?tid=1206743216#/learn/content?type=detail&id=1211570840&cid=1214317955

一、数值微分与数值积分

数值微积分不同于高等数学中的微积分，高等数学是对函数采用解析方法求解，数值微积分适合求解没有或很难求出微分或积分表达式的问题的计算。

1.1 数值微分

当h非常小时，得到：

函数f(x)在$x_0$点的微分，接近于函数在该点的差分，而f(x)在$x_0$的导数，接近于函数在该点的差商。

%% 6.1 数值微分与数值积分
%例1  设f(x)=sinx，在[0，2π]范围内随机采样，计算f'(x)的近似值，并与理论值f'(x)=cosx 进行比较。
x=[0,sort(2*pi*rand(1,5000)),2*pi];
y=sin(x);
f1=diff(y)./diff(x);
f2=cos(x(1:end-1));
plot(x(1:end-1),f1,x(1:end-1),f2);
d=norm(f1-f2)  %范数

1.2 数值积分

当被积函数的原函数无法用初等函数表示时，或被积函数是用表格形式给出的，这时就需要用数值解法来求定积分的近似值，求定积分的数值方法多种多样，如梯形法、辛普森法、高斯求积公式等。

在每一个小的区间上，定积分的值可以近似求得，从而避免了牛顿莱布尼兹公式需要寻求原函数的困难。

quad函数

quadl函数

%例2  分别用quad函数和quadl函数求定积分的近似值，并在相同的积分精度下，比较被积函数的调用次数。
format long
f=@(x) 4./(1+x.^2);
[I,n]=quad(f,0,1,1e-8)
 
[I,n]=quadl(f,0,1,1e-8)
 
(atan(1)-atan(0))*4
 
format short

integral函数

%例3 求定积分。 
%被积函数文件fe.m：
function f=fe(x)
f=1./(x.*sqrt(1-log(x).^2));

I=integral(@fe,1,exp(1))

quadgk函数

求振荡函数的定积分

%例4  求定积分。
%被积函数文件fsx.m：
function f=fsx(x)
f=sin(1./x)./x.^2;

I=quadgk(@fsx,2/pi,+Inf)

**当函数关系表达式不知道，只有实验测定的一组样本点和样本值，这样就无法使用quad等函数来计算定积分。**针对这种情况，有以下方法进行求解：

trapz函数

%例5  设x=1:6，y=[6,8,11,7,5,2],用trapz函数计算定积分。
x=1:6;
y=[6,8,11,7,5,2];
plot(x,y,'-ko');
grid on
axis([1,6,0,11]);
I1=trapz(x,y)

I2=sum(diff(x).*(y(1:end-1)+y(2:end))/2)   %梯形面积公式

多重定积分的数值求解

% 例6  分别求二重积分和三重积分。 
f1=@(x,y) exp(-x.^2/2).*sin(x.^2+y);
I1=quad2d(f1,-2,2,-1,1)

f2=@(x,y,z) 4*x.*z.*exp(-z.*z.*y-x.*x);
I2=integral3(f2,0,pi,0,pi,0,1)

二、线性方程组求解

2.1 直接法

就是在不考虑舍入误差的情况下，通过有限步的矩阵初等运算，来求得方程组精确解的方法。

2.1.1 高斯（Gauss）消去法

2.1.2 列主元消去法

是目前计算机上求解线性方程组的标准算法，其特点就是通过消元，将一般线性方程组的求解问题转化为三角方程组的求解问题。

A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';  %b是列向量
x=A\b

可以用矩阵求逆函数来求解线性方程组，即x等于A的逆矩阵乘以b,但在实际计算中，不建议使用求逆矩阵的方法，因为求逆方法需要更多的计算量，它先求逆矩阵，在求矩阵乘法，而左除运算就只有一次除法，况且，求逆算法给出的可能还是一个不精确的答案。

2.1.3 矩阵的三角分解法

矩阵分解是设计算法的重要技巧，是指将一个给定的矩阵分解成若干个特殊类型矩阵的乘积，从而将一个一般的矩阵计算问题，转化为几个易求的特殊矩阵的计算问题。这种方法的优点是运算速度快，可以节省存储空间。

LU分解

将一个n阶矩阵表示为一个下三角矩阵和一个上三角矩阵的乘积，线性代数中已证明，只要方阵是非奇异的，LU分解总是可以进行的。

对于三角矩阵很容易求解，可以先求解向量y使Ly等于b，再求解Ux等于y

使用第一种方式，矩阵L往往不是下三角阵，但可以通过行交换，成为一个下三角阵。

%例2  用LU分解求解例1中的线性方程组。
A=[2,1,-5,1;1,-5,0,7;0,2,1,-1;1,6,-1,-4];
b=[13,-9,6,0]';
[L,U]=lu(A);
x=U\(L\b)
%第二种格式
[L,U,P] = lu(A);
x = U\(L\P*b)

QR分解

Cholesky分解(乔累斯分解)

2.2 迭代法

就是先给定一个解的初始值，然后按照一定的迭代算法进行逐步逼近，求出更精确的近似解的方法。
不断用变量的原值推出它的新值的过程，是用计算机解决问题的一种基本方法。

2.2.1 雅克比（Jacobi）迭代法

对于线性方程组Ax=b，如果A为非奇异方阵且主对角线元素不等于0，则可将A分解为D-L-U，其中D为对角阵，其元素为A的对角元素，L与U为A的下三角阵取反和上三角阵取反。

%雅可比迭代法的函数文件jacobi.m：
function [y,n]=jacobi(A,b,x0,ep)
D=diag(diag(A));
L=-tril(A,-1);
U=-triu(A,1