数值微分与数值积分
1. 数值微分
1) 求向前差分的函数diff( )
调用格式
- 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)=diff(diff(x)) - dx=diff(A,n,dim):计算矩阵A的n阶差分,dim=1时(默认状态),按列计算差分;dim=2时,按行计算差分。
e.g.设f(x)=sinx,在[0,2pi]范围内随机采样,计算f’(x)的近似值,并与理论值f’(x)=cos(x)进行比较。
>> x=[0,sort(2*pi*rand(1,5000)),2*pi];
>> y=sin(x);
>> f1=diff(y)./diff(x);
%f'(x)的表示方法
>> f2=cos(x(1:end-1));
>> plot(x(1:end-1),f1,x(1:end-1),f2);
>> d=norm(f1-f2)
%范数norm(A)返回A的最大奇异值,即max(svd(A))
d =
0.0456
2. 数值积分
1)数值积分基本原理
2)数值积分的实现
i. 基于全局自适应积分方法
I=integral(filename,a,b)
I:定积分值
filename:被积函数
a,b:定积分下限、上限,积分限可以无穷大
e.g.求定积分
%先定义被积函数文件fe.m
function f = fe(x)
f = 1./(x.*sqrt(1-log(x).^2));
%注意使用点除,点乘方
end
%命令行
>> I=integral(@fe,1,exp(1)) %文件名
I =
1.5708
ii. 基于自适应辛普森方法
[I,n]=quad(filename,a,b,tol,trace)
iii. 基于自适应Gauss-Lobatto方法
[I,n]=quadl(filename,a,b,tol,trace)
filename:被积函数名;
a,b为定积分的下限和上限,积分限[a,b]必须是有限的,不能为无穷大(inf)
tol:控制积分精度,默认取 tol=10(-6)
trace:控制是否展现积分过程,若取非0则展现积分过程,取0则不展现,默认取0
I:返回参数,即定积分值
n:被积函数的调用次数
e.g.分别用quad函数和quadl函数求定积分的近似值,并在相同的积分精度下,比较被积函数的调用次数。
>> format long
>> f=@(x) 4./(1+x.^2);
%定义函数方式
>> [I,n]=quad(f,0,1,1e-8)
I =
3.141592653733437
n =
61
>> [I,n]=quadl(f,0,1,1e-8)
I =
3.141592653589806
n =
48
>> format short
iv. 基于自适应高斯-克朗罗德方法
[I,err]=quadgk(filename,a,b)
err:近似误差范围
积分上下限可以是无穷大(-Inf或Inf),也可以是复数
若积分上下限为复数,则quadgk( )函数在复平面上求积分
e.g.求定积分
被积函数文件fe.m
function f = fe(x)
f = sin(1./x)./(x.^2);
end
命令行
>> I=quadgk(@fe,2/pi,+Inf) %文件名
I =
1.0000
v. 基于梯形积分法
已知
(
x
i
,
y
i
)
(x_i,y_i)
(xi,yi)(i=1,2,…,n),且
a
=
x
1
<
x
2
<
…
<
x
n
=
b
a=x_1<x_2<…<x_n=b
a=x1<x2<…<xn=b,求
I
=
∫
a
b
f
(
x
)
d
(
x
)
I=∫_a^bf(x)d(x)
I=∫abf(x)d(x)近似值。I=trapz(x,y),其中,向量x,y定义函数关系y=f(x)。
原理:trape( )函数采用梯形积分法则,积分近似值为:
可用以下语句实现:sum(diff(x).*(y(1:end-1)+y(2:end))/2)
e.g.设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]);
%设值x轴数值范围为1-6,y轴数值范围为0-11
>> I1=trapz(x,y)
I1 =
35
>> I2=sum(diff(x).*(y(1:end-1)+y(2:end))./2)
I2 =
35
3. 多重定积分的数值求解
1)求二重积分的数值解 ∫ c d ∫ a b f ( x , y ) d x d y ∫_c^d∫_a^bf(x,y)dxdy ∫cd∫abf(x,y)dxdy
I=integral2(filename,a,b,c,d)
I=quad2d(filename,a,b,c,d)
I=dblquad(filename,a,b,c,d,tol)
2)求三重积分的数值解 ∫ e f ∫ c d ∫ a b f ( x , y , z ) d x d y d z ∫_e^f∫_c^d∫_a^bf(x,y,z)dxdydz ∫ef∫cd∫abf(x,y,z)dxdydz
I=integral3(filename,a,b,c,d,e,f)
I=triplequad(filename,a,b,c,d,e,f,tol)
e.g.分别求二重积分和三重积分
>> f1=@(x,y) exp(-x.^2/2).*sin(x.^2+y);
>> I1=quad2d(f1,-2,2,-1,1)
I1 =
1.5745
>> f2=@(x,y,z) 4*x.*z.*exp(-z.^2.*y-x.^2);
>> I2=integral3(f2,0,pi,0,pi,0,1)
I2 =
1.7328
%注意:1)函数定义方式;2)函数参数个数;3)积分下限、上限书写顺序