第六节-数值微分与方程求解

01-数值微分与数值积分




% 第六节-数值微分与方程求解
% 01-数值微分与数值积分


% 数值微分 和解析几何 不同

%(1)数值差分与差商
%   向前差分:△f(xo)= f(xo+ h)- f(xo)
%   向后差分: ▽f(xo)= f(xo)- f(xo-h)
%   中心差分: δf(x0)= f(xo+ h/2)- f(xo- h/2)

%   向前差商:f'(x0)=△f(x0)/h
%   向后差商:f'(x0)=▽f(x0)/h
%   中心差商:f'(x0)=δf(x0)/h

%   MATLAB提供了求向前差分的函数diff,其调用格式有三E种: 
%   dx=diff(x): 计算向量x的一-阶向前差分，dx(i)=x(i+1)-x(i)，i=1,2，...，n-1。
%   dx=diff(x,n): 计算向量x的n阶向前差分。
%   例如，diff(x,2)=dif(diff(x))。
%   dx=difflA,n,dim): 计算矩阵A的n阶差分，dim=1时 (默认状态)
%   按列计算差分; dim=2， 按行计算差分。


%   例1设f(x)=sin x,在[0， 2 π ]范围内随机采样，计算f(x)的近似值，并与
%   理论值f(x)=cos x进行比较。
%   
x=[0,sort(2*pi*rand(1 ,5000)),2*pi]; %再【0；2Π】内进行随机采样
% 对每一个数据进行随机采样，一共5000个，随后对这个随机数进行由小到大排列
y=sin(x);
f1=diff(y)./ diff(x); %求一阶向前差分，得到差商向量
f2=cos(x(1:end-1));   %求各地导数的理论值f2
plot(x(1:end-1),f1 ,x(1:end-1),f2); %绘制近似值和理论值直线，发现两者较为接近
d=norm(f1-f2)   %求近似值和理论值之间的范数


% (1)数值积分基本原理
% 牛顿一莱布尼兹( Newton-Leibniz)公式:
% 数值积分
% 将积分区间[a，b]分成n个子区间[xi，xi+1], i=1, 2，...，n，其中
% x=a，Xn+1=b, 这样求定积分问题就分解为下面的求和问题:
% S=∑f(x) 积分变成求和

%数值积分
% (2)数值积分的实现
%   基于自适应辛普森方法
%   [l,n]=quadl(filename,a,b,tol,trace)
%   基于自适应Gauss-Lobatto方法
%   [I,n]=quadl(filename,a,b,tol,trace)

% 其中，filename是被积函数名; 
% a和b分别是定积分的下限和上限，
% 积分限[a，b]必须是有限的，不能为无穷大(Inf) ;
% tol用来控制积分精度，默认时取tol=10-6; 
% trace控制是否展现积分过程，若取非0则展现积分过程
% 取0则不展现，默认时取trace=0;返回参数即定积分的值，
% n为被积函数的调用次数。

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

%数值积分的实现
%   基于全局自适应积分方法
%   l=integral(filename,a,b)

% (2)数值积分的实现
% 基于自适应高斯-克朗罗德方法
% [I,err]=quadgkfilename,a,b)
% 其中，err返回近似误差范围，其他参数的含义和用法与quad
% 函数相同。积分.上下限可以是无穷大( -Inf或Inf)，也可以是复数。
% 如果积分上下限是复数，则quadgk函数在复平面.上求积分。

%数值积分的实现
%基于梯形积分法
%已知(x，y) (i=1, 2, ..，n)，且a=x<x*<...<Xn=b,求1=[” f(x)dx近似值。
% I=trapz(x,y)


% eg 
x1=[1:6]
y1=[6,8,11,7,5,2]
figure(2)
plot(x1,y1,'-ko')
grid on
axis([1,6,0,11])
format 
I1=trapz(x1,y1)  %使用梯形计算法求积分
I2=sum(diff(x1).*(y1(1:end-1)+y1(2:end))/2) %q求和公式求出各个梯形面积，最后求和
% I1=I2


%多重定积分数值求解
%求二重积分的数值解:求积分f(x,y)dxdy
%   l=integral2(filename,a,b,c,d)
%   l=quad2d(filename,a,b,c,d)
%   l=dblquad(filename,a,b,c,d,tol)
% 求三重积分的数值解:S (x,y,z)dxdydz
%   l=integral3(filename,a,b,c,d,e,f)
%   l=triplequad(filename,a,b,c,d,e,f,tol)



% eg 求y=[-1,1] x=[-2,2] e^(-x^2/2)sin(x^2 +y)dxdy 的二重积分
f1=@(x,y)exp(-x.^2/2).*sin(x.^2+y)
I1=quad2d(f1,-2,2,-1,1)

% eg 求 x=[0:Π] y=[0:Π] z=[0:1] 4xze^(-z^2y-x^2)dxdydz 的三重积分
f2=@(x,y,z)4*x.*z.*exp(-z.*z.*y-x.*x)
I2=integral3(f2,0,pi,0,pi,0,1)

更多文章请关注《万象专栏》