备战数学建模7-MATLAB数值微积分与方程求解
Posted nuist__NJUPT
tags:
篇首语:本文由小常识网(cha138.com)小编为大家整理,主要介绍了备战数学建模7-MATLAB数值微积分与方程求解相关的知识,希望对你有一定的参考价值。
目录
一、数值微分与数值积分
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));
绘制的图形如下所示:
以上是关于备战数学建模7-MATLAB数值微积分与方程求解的主要内容,如果未能解决你的问题,请参考以下文章