备战数学建模7-MATLAB数值微积分与方程求解

一、数值微分与数值积分

1-数值微分的基本原理

2-数值微分的实现

3-数值积分的基本原理

4-数值积分的实现方法

二、线性方程组

三、非线性方程组

1-单变量非线性方程的求解

2-函数极值的计算

3-最小值问题的应用实例

四、常微分方程

1-常微分方程的一般概念

2-常微分方程的求解函数

一、数值微分与数值积分

1-数值微分的基本原理

当f(x)函数比较复杂，或者以离散列表的形式给出的时候，如果要求导数值，可以用数值方法。

任意函数f(x)在x0点的导数是通过极限定义的，有如下三种形式：

如果去掉极限h趋向于0的极限过程，我们可以理解为该函数是以h为步长的向前，向后，中心差分。差分近似等x在x0点的微分，如下所示：

下面我们看一下差商公式，近似等于x在x0点的的导数。

2-数值微分的实现

MATLAB中提供了求向前差分的函数diff()，其有如下三种调用格式：

(1) dx = diff(x) 表示计算向量x的一阶向前差分，dx(i) = x(i+1) - x(i).

(2) dx = diff(x,n) 表示计算x的n阶向前差分，diff(x,2) = diff(diff(x)).

(3) dx = diff(x,n.,dim) 表示计算矩阵A的n阶差分，dim=1,是默认差分，案列计算差分，dim=2时，按行计算差分.

我们看一下上面的例子1，具体的代码如下：

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) ; %绘制f1和f2的曲线
d = norm(f1 - f2) %求范数

或者的图像我们可以看出，两条曲线基本重合，所以f1和f2近似相等。

3-数值积分的基本原理

如果函数的原函数可以找到，那么可以使用牛顿-莱布尼茨公式求解定积分，如下所示：

但是有些时候原函数很难找到，这时候就需要使用数值积分的方式求解定积分。

数值积分的基本原理就是将区间划分为许多的小区间，这样求定积分就分解成求和问题，如下所示：

4-数值积分的实现方法

MATLAB中提供了求定积分的方法，函数的调用格式如下：

1-基于自适应的辛普森方法，[I,n] = quad(filname,a,b,tol,trace)

2-基于自适应的Gauss-Labatto方法 [I,n] = quadl(filename,a,b,tol,trace)

其中，filename是被积函数名，a,b是积分上下限，tol用来控制积分精度，trace控制是否展现积分过程，返回的I是定积分的值，返回的n是被调用的次数。

我们看一下上面的例子2，我们通过调用两个方法计算的结果I都是等于3.14，而quad()函数的调用次数是61次，而quadl()函数的调用次数是48次，故此处明显 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)

MATLAB中也提供了如下求解定积分的方法，基于全局自适应的球定积分方法，

I = integral(filename, a, b) ,其中，filename是被积函数，a,b是定积分的上下限，积分限可以无穷大。

我们看一下上面的例子3，求解定积分。

f = @(x) 1./(x.*sqrt(1-log(x).^2)) ;
I = integral(f,1,exp(1))

MATLAB中还提供了基于高斯-克朗罗德方法去求解震荡函数的定积分，函数的调用格式如下：

[I,err] = quadgk(filename, a, b) ，其中err返回近似误差范围，积分上下限可以是无穷大，也可以是复数，即在复数平面去积分

我们看一下上面的求震荡函数定积分的例子4，具体的代码如下所示:

f = @(x) sin(1./x)./(x.^2);
I = quadgk(f, 2/pi, +inf)

MATLAB中提供对于多重积分的求解方法，具体如下所示：

我们看一下上面的例子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)

二、线性方程组

一般包含两大类方法，即直接求解发和迭代法。

我们先看第一类，直接求解法：

1-高斯消去法；2-列主元消去法；3-矩阵的三角分解法

利用左除运算符直接求解，如下所示：

我们看一下上面的例子1，用直接法求解线性方程组，具体代码如下：

A = [2,1,-5,1; 1,-5,0,7; 0,2,1,-1; 1,6,-1,-4] ;
b = [13; -9; 6; 0] ;
A\b

下面看一下利用矩阵分解方法求解线性方程组。

MATLAB的LU分解函数的调用格式如下：

[L,U] = lu(A)表示产生一个上三角矩阵U和一个变换形式的下三角矩阵L，使得满足A=LU，注意，这里的矩阵A必须是方阵。

[L,U,P] = lu(A) 表示产生一个上三角矩阵L和一个下三角矩阵U和一个置换矩阵P，使得PA=LU

具体的转换形式如下所示，两种矩阵分解方法求线性方程组。

我们看一下上面的例子2，使用LU分解方法求解线性方程组，代码如下：

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)

三、非线性方程组

1-单变量非线性方程的求解

MATLAB中提供fzero()函数求解单变量非线性方程，函数的调用格式入下所示：

fzero(filename, x0) ,其中filename是方程左端的函数表达式，x0是初始值。

我们看一下上面的例子1，用fzero()函数求解非线性方程组，迭代初始值的选取很重要。

代码如下：

f = @(x) x - 1./x + 5 ;
x1 = fzero(f, -5)
x2 = fzero(f, 1)

MATLAB中提供了fsolve()函数完成非线性方程组的求解，函数的具体调用个数如下：

x = fsolve(filename, x0, option),x为返回的近似解，filename为待求根方程左端的表达式，x0是初值，option用于设置优化工具箱的参数，可以调用optimset函数来设置。

我们可以看一下上面的例子2，使用fsolve()函数求解非线性方程组在(1,1,1)附近的解，代码如下：

 f = @(x) [sin(x(1)) + x(2)+x(3)^2*exp(x(1)), x(1)+x(2)+x(3),x(1)*x(2)*x(3)] ;x = fsolve(f, [1,1,1], optimset('Display', 'off'))

2-函数极值的计算

MATLAB中秩考虑函数极小值问题，如果要求极大值，可以用求函数负值的极小值来解决。

1-无约束优化

MATLAB提供了三种求极小值的函数。

[xmin, fmin] = fminbnd(filename, x1, x2, option)

[xmin, fmin] = fminsearch(filename, x0,option)

[xmin, fmin] = fminunc(filename,x0,option)

其中filename是定义的目标函数，第一个函数的x1，x2表示被研究的区间的左右边界，后面两个函数的x0是一个向量，表示极值点的初值，option表示优化选项。

我们看一下上面的例子3，求函数在求间内的最小值点，代码如下所示：

 f = @(x) x - 1./x + 5 ;[xmin, fmin] = fminbnd(f, -10, -1) [xmin, fmin] = fminbnd(f, 1, 10)

2-有约束的优化

MATLAB中提供了求解有约束的条件下的函数最小值，函数的调用格式为：

[xmin, fmin] = fmincon(filename, x0,A,b,Aeq,beq,Lbnd,Ubnd,Nonf,option)

其中filename为函数，x0为初始值，A为系数矩阵，b为结果矩阵，其余分别为线性不等式约束，线性等式约束，x的下界和上界，非线性约束，优化选项。如果约束条件不存在，可以用空矩阵来表示。

我们使用fmincon()函数求解有约束的最小值问题。

f =@(x) 0.4*x(2) + x(1)^2 + x(2)^2 - x(1)*x(2)  + 1/30*x(1)^3 ;
A = [-1,-0.5; -0.5,-1] ;
b = [-0.4; -0.5] ;
lb = [0;0] ;
x0 = [0.5; 0.5] ;
option = optimset('Display','off') ;
[xmin, fmin] = fmincon(f, x0, A,b, [], [], lb, [], [], option)

3-最小值问题的应用实例

我们看一下上面的例子5，其实就是求方程式最小值问题，由题意知：

假设仓库所选点的坐标为(x,y),则总里程的表达式为：

代码如下所示，求出在坐标点(19.814348806223666 41.124724254226052)处建立建立仓库，每周货车的总里程最少，最少为 1361.8公里。

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])

我们看一下例子6，这个是加了一个非线性约束条件，我们可以建立一个约束函数funny,代码如下：

function [c,h] = funny(x)
c = [] ;
h = x(2) - x(1)^2 ;
end

在命令行输入如下代码，完成计算，代码如下：

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] = fmincon(f,[15,30],[],[],[],[],[],[],'funny')

四、常微分方程

1-常微分方程的一般概念

求常微分方程初值问题，就是寻找函数y(t)，使其满足如下方程：

2-常微分方程的求解函数

MATLAB中给出了求解常微分方程的函数，其函数调用格式入下：

[t,y] = solver(filename, tspan, y0, option), 其中，t和y分别给出时间向量和数值解，filename是定义f(t,y)的函数名，该函数必须返回一个列向量，tspan表示求解区间，y0是初始向量，option是可选参数，用于设置求解属性。

另外常微分方程的数值求解函数的统一命名格式如下所示：

我们看一下上面的例子1，代码如下所示，并将近似值和精确值绘成曲线继续比较。

f = @(t,y) (y^2-t-2)/(4*(t+1)) ;
[t,y] = ode23(f,[0,10],2) ;
y1 = sqrt(t+1) + 1 ;
plot(t,y,'b:',t,y1,'r');

绘制的图形如下所示：

对于高阶常微分方程，需要先转换为一阶常微分方程，然后进行求解。

对于上面的例子2，改成一阶的常微分方程后，代码如下所示：

f = @(t,x) [-2,0; 0,1] *[x(2); x(1)] ;
[t,x] = ode45(f,[0,20],[1,0]) ;
subplot(2,2,1) ;
plot(t,x(:,2)) ;
subplot(2,2,2)
plot(x(:,2), x(:,1));

绘制的图形如下所示：