Matlab与工程计算 第七章 数值积分、微分
合集下载
相关主题
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
oad wind div = divergence(x,y,z,u,v,w); slice(x,y,z,div,[90 134],[59],[0]); shading interp daspect([1 1 1]) camlight
Xiamen University
Matlab and Engineering Calculation
dx=diff(f([x,3.01]))/0.01; %直接对f(x)求数值导数 gx=g(x); %求函数f的导函数g在假设点的导数
plot(x,dpx,x,dx,'.',x,gx, 'r-'); %作图
Xiamen University
Matlab and Engineering Calculation
f ' ( x) =
f ( x + h) − f ( x ) h
Xiamen University
Matlab and Engineering Calculation
例6-6 生成以向量V=[1,2,3,4,5,6]为基础的范得蒙矩阵，按列 进行差分运算。 命令如下： V=vander(1:6) DV=diff(V) %计算V的一阶差分
7.1.2 数值积分的实现方法(Quadrature)
1．自适应辛普生法 基于变步长辛普生法，MATLAB给出了quad函数来 求定积分。该函数的调用格式为： [I,n]=quad(fun,a,b,tol,trace) fun是被积函数名、句柄或内联函数对象; a和b分别是定积分的下限和上限; tol用来控制积分精度，缺省时取tol=10-6; trace控制是否展现积分过程，若取非0则展现积分过 程，取0则不展现，缺省时取trace=0 I即定积分值，n为被积函数的调用次数。
标量场梯度计算(gradient)
[Fx,Fy,Fz,...] = gradient(F,h1,h2,...) [Fx,Fy,Fz,...] = gradient(F)
v = -2:0.2:2; [x,y] = meshgrid(v); z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.2); contour(v,v,z), hold on, quiver(v,v,px,py), hold off
Xiamen University
Xiamen University
Matlab and Engineering Calculation
例6-7 用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出 f'(x)的图像。 程序如下： f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2'); g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5'); x=-3:0.01:3; p=polyfit(x,f(x),5); dp=polyder(p); dpx=polyval(dp,x); %用5次多项式p拟合f(x) %对拟合多项式p求导数dp %求dp在假设点的函数值
调用函数quadl求定积分： [I,n]=quadl(fx,1,2.5,1e-10) I= 0.28579444254754 n= 33
Xiamen University
Matlab and Engineering Calculation
3．被积函数由一个表格定义
在MATLAB中，对由表格形式定义的函数关 系的求定积分问题用trapz(X,Y)函数。其中向 量X,Y定义函数关系Y=f(X)。
Xiamen University
Matlab and Engineering Calculation
例7-2 求定积分。 (1) 被积函数文件fx.m。 f=inline('x.*sin(x)./(1+cos(x).*cos(x)) '); I=quad8(f,0,pi) I= 2.4674
Xiamen University
Xiamen University
Matlab and Engineering Calculauad8和quadl函数求定积分的近似 值，并在相同的积分精度下，比较函数的调用次数。
调用函数quad求定积分： format long; fx=inline('exp(-x)'); [I,n]=quad(fx,1,2.5,1e-10) I= 0.28579444254766 n= 65 调用函数quad8求定积分： [I,n]=quad8(fx,1,2.5,1e-10) I= 0.28579444254754 n= 33
Xiamen University
Matlab and Engineering Calculation
2．牛顿－柯特斯法
基于牛顿－柯特斯法，MATLAB给出了 quad8函数来求定积分。该函数的调用格式为： [I,n]=quad8(fun,a,b,tol,trace) 其中参数的含义和quad函数相似。该函数可 以更精确地求出定积分的值，且一般情况下函数 调用的步数明显小于quad函数，从而保证能以更 高的效率求出所需的定积分值。
Matlab and Engineering Calculation
拉普拉斯算子Discrete Laplacian
L = del2(F) L = del2(F,h) L = del2(F,hx,hy) L = del2(F,hx,hy,hz,...)
[x,y] = meshgrid(-4:4,-3:3); U = x.*x+y.*y; V = 4*del2(U);
Xiamen University
Matlab and Engineering Calculation
7.2 数值微分
7.2.1 数值差分与差商 7.2.2 数值微分的实现 在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函 数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，按行计算差分。
Xiamen University
Matlab and Engineering Calculation
例7-1 求定积分。 (1) 建立被积函数文件fesin.m。 function f=fesin(x) f=exp(-0.5*x).*sin(x+pi/6); (2) 调用数值积分函数quad求定积分。 [S,n]=quad('fesin',0,3*pi) S= 0.9008 n= 77
7.1.3 二重定积分的数值求解 使用MATLAB提供的dblquad函数就可 以直接求出上述二重定积分的数值 解。该函数的调用格式为： I=dblquad(f,a,b,c,d,tol,trace) 该函数求f(x,y)在[a,b]×[c,d]区域上的二 重定积分。参数tol，trace的用法与函 数quad完全相同。
Xiamen University
Matlab and Engineering Calculation
已知一维离散函数，求其导数 x=0:0.1:10; y=sin(x); dydx1=gradient(y,0.1); dydx2=diff(y)/0.1;
⎧ f ( x + h) − f ( x ) ⎪ ⎪ h ' f ( x) = ⎨ ⎪ f ( x + h) − f ( x − h) ⎪ 2h ⎩
Xiamen University
Matlab and Engineering Calculation
7.1.2 数值积分的实现方法(Quadrature)
3
2.5
2
f(x)
1.5
1
0.5
0
0
1
2
3 x
4
5
6
Xiamen University
Matlab and Engineering Calculation
plot(x,y,x,dydx1,x(1:end-1),dydx2, '+')
Xiamen University
Matlab and Engineering Calculation
矢量场散度计算
div = divergence(X,Y,Z,Fx,Fy,Fz)
∂Fx ∂Fy ∂Fz ∇⋅F = + + ∂x ∂y ∂z
Matlab and Engineering Calculation
第7章 MATLAB数值积分与微分 7.1 数值积分 7.2 数值微分
Xiamen University
Matlab and Engineering Calculation
7.1 数值积分
7.1.1 数值积分基本原理 求解定积分的数值方法多种多样，如简单的 梯形法、辛普生(Simpson)法、牛顿－柯特斯 (Newton-Cotes)法等都是经常采用的方法。 它们的基本思想都是将整个积分区间[a,b]分 成n个子区间[xi,xi+1]，i=1,2,…,n，其中x1=a， xn+1=b。这样求定积分问题就分解为求和问题。
b ⎧ I = ∫ f ( x ' )dx ' ⎪ a ⎪ x ⎪ ' ' ⎨ I ( x) − I (a ) = ∫a f ( x )dx ⎪ ⎪ f ( x) = dI ( x) ⎪ dx ⎩
Xiamen University
Matlab and Engineering Calculation
Matlab and Engineering Calculation
2．高阶方法
MATLAB给出了quadl函数来替代quad8求定积分 [I,n]=quadl(fun,a,b,tol,trace) 其中参数的含义和quad函数相似。该函数可 以更精确地求出定积分的值，且一般情况下函数 调用的步数明显小于quad函数，从而保证能以更 高的效率求出所需的定积分值。
Xiamen University
Matlab and Engineering Calculation
例7-5 计算二重定积分 (1) 建立一个函数文件fxy.m： function f=fxy(x,y) global ki; ki=ki+1; %ki用于统计被积函数的调用次数 f=exp(-x.^2/2).*sin(x.^2+y); (2) 调用dblquad函数求解。 global ki;ki=0; I=dblquad('fxy',-2,2,-1,1) ki I= 1.57449318974494 ki = 1038
Trapezoidal numerical integration Z = trapz(X,Y) Z = cumtrapz(X,Y)
例7-4 已知函数在离散点的值，计算定积分值和积 分函数。 命令如下： x=pi/4:0.01:2*pi; y=sin(-x); I=trapz(x,y) Ix=cumtrapz(x,y) plot(x,y,’-o’,x,Ix)
矢量场旋度计算
[curlx,curly,curlz] = curl(X,Y,Z,U,V,W)
⎛ ∂Fz ∂Fy ⎞ ∇× F = i⎜ − ⎟+ ∂z ⎠ ⎝ ∂y ⎛ ∂Fx ∂Fz j⎜ − ∂x ⎝ ∂z ⎛ ∂Fy ∂Fx ⎞ ⎞ − ⎟ ⎟+k⎜ ∂x ∂y ⎠ ⎠ ⎝
Xiamen University
Xiamen University
Matlab and Engineering Calculation
dx=diff(f([x,3.01]))/0.01; %直接对f(x)求数值导数 gx=g(x); %求函数f的导函数g在假设点的导数
plot(x,dpx,x,dx,'.',x,gx, 'r-'); %作图
Xiamen University
Matlab and Engineering Calculation
f ' ( x) =
f ( x + h) − f ( x ) h
Xiamen University
Matlab and Engineering Calculation
例6-6 生成以向量V=[1,2,3,4,5,6]为基础的范得蒙矩阵，按列 进行差分运算。 命令如下： V=vander(1:6) DV=diff(V) %计算V的一阶差分
7.1.2 数值积分的实现方法(Quadrature)
1．自适应辛普生法 基于变步长辛普生法，MATLAB给出了quad函数来 求定积分。该函数的调用格式为： [I,n]=quad(fun,a,b,tol,trace) fun是被积函数名、句柄或内联函数对象; a和b分别是定积分的下限和上限; tol用来控制积分精度，缺省时取tol=10-6; trace控制是否展现积分过程，若取非0则展现积分过 程，取0则不展现，缺省时取trace=0 I即定积分值，n为被积函数的调用次数。
标量场梯度计算(gradient)
[Fx,Fy,Fz,...] = gradient(F,h1,h2,...) [Fx,Fy,Fz,...] = gradient(F)
v = -2:0.2:2; [x,y] = meshgrid(v); z = x .* exp(-x.^2 - y.^2); [px,py] = gradient(z,.2,.2); contour(v,v,z), hold on, quiver(v,v,px,py), hold off
Xiamen University
Xiamen University
Matlab and Engineering Calculation
例6-7 用不同的方法求函数f(x)的数值导数，并在同一个坐标系中做出 f'(x)的图像。 程序如下： f=inline('sqrt(x.^3+2*x.^2-x+12)+(x+5).^(1/6)+5*x+2'); g=inline('(3*x.^2+4*x-1)./sqrt(x.^3+2*x.^2-x+12)/2+1/6./(x+5).^(5/6)+5'); x=-3:0.01:3; p=polyfit(x,f(x),5); dp=polyder(p); dpx=polyval(dp,x); %用5次多项式p拟合f(x) %对拟合多项式p求导数dp %求dp在假设点的函数值
调用函数quadl求定积分： [I,n]=quadl(fx,1,2.5,1e-10) I= 0.28579444254754 n= 33
Xiamen University
Matlab and Engineering Calculation
3．被积函数由一个表格定义
在MATLAB中，对由表格形式定义的函数关 系的求定积分问题用trapz(X,Y)函数。其中向 量X,Y定义函数关系Y=f(X)。
Xiamen University
Matlab and Engineering Calculation
例7-2 求定积分。 (1) 被积函数文件fx.m。 f=inline('x.*sin(x)./(1+cos(x).*cos(x)) '); I=quad8(f,0,pi) I= 2.4674
Xiamen University
Xiamen University
Matlab and Engineering Calculauad8和quadl函数求定积分的近似 值，并在相同的积分精度下，比较函数的调用次数。
调用函数quad求定积分： format long; fx=inline('exp(-x)'); [I,n]=quad(fx,1,2.5,1e-10) I= 0.28579444254766 n= 65 调用函数quad8求定积分： [I,n]=quad8(fx,1,2.5,1e-10) I= 0.28579444254754 n= 33
Xiamen University
Matlab and Engineering Calculation
2．牛顿－柯特斯法
基于牛顿－柯特斯法，MATLAB给出了 quad8函数来求定积分。该函数的调用格式为： [I,n]=quad8(fun,a,b,tol,trace) 其中参数的含义和quad函数相似。该函数可 以更精确地求出定积分的值，且一般情况下函数 调用的步数明显小于quad函数，从而保证能以更 高的效率求出所需的定积分值。
Matlab and Engineering Calculation
拉普拉斯算子Discrete Laplacian
L = del2(F) L = del2(F,h) L = del2(F,hx,hy) L = del2(F,hx,hy,hz,...)
[x,y] = meshgrid(-4:4,-3:3); U = x.*x+y.*y; V = 4*del2(U);
Xiamen University
Matlab and Engineering Calculation
7.2 数值微分
7.2.1 数值差分与差商 7.2.2 数值微分的实现 在MATLAB中，没有直接提供求数值导数的函数，只有计算向前差分的函 数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，按行计算差分。
Xiamen University
Matlab and Engineering Calculation
例7-1 求定积分。 (1) 建立被积函数文件fesin.m。 function f=fesin(x) f=exp(-0.5*x).*sin(x+pi/6); (2) 调用数值积分函数quad求定积分。 [S,n]=quad('fesin',0,3*pi) S= 0.9008 n= 77
7.1.3 二重定积分的数值求解 使用MATLAB提供的dblquad函数就可 以直接求出上述二重定积分的数值 解。该函数的调用格式为： I=dblquad(f,a,b,c,d,tol,trace) 该函数求f(x,y)在[a,b]×[c,d]区域上的二 重定积分。参数tol，trace的用法与函 数quad完全相同。
Xiamen University
Matlab and Engineering Calculation
已知一维离散函数，求其导数 x=0:0.1:10; y=sin(x); dydx1=gradient(y,0.1); dydx2=diff(y)/0.1;
⎧ f ( x + h) − f ( x ) ⎪ ⎪ h ' f ( x) = ⎨ ⎪ f ( x + h) − f ( x − h) ⎪ 2h ⎩
Xiamen University
Matlab and Engineering Calculation
7.1.2 数值积分的实现方法(Quadrature)
3
2.5
2
f(x)
1.5
1
0.5
0
0
1
2
3 x
4
5
6
Xiamen University
Matlab and Engineering Calculation
plot(x,y,x,dydx1,x(1:end-1),dydx2, '+')
Xiamen University
Matlab and Engineering Calculation
矢量场散度计算
div = divergence(X,Y,Z,Fx,Fy,Fz)
∂Fx ∂Fy ∂Fz ∇⋅F = + + ∂x ∂y ∂z
Matlab and Engineering Calculation
第7章 MATLAB数值积分与微分 7.1 数值积分 7.2 数值微分
Xiamen University
Matlab and Engineering Calculation
7.1 数值积分
7.1.1 数值积分基本原理 求解定积分的数值方法多种多样，如简单的 梯形法、辛普生(Simpson)法、牛顿－柯特斯 (Newton-Cotes)法等都是经常采用的方法。 它们的基本思想都是将整个积分区间[a,b]分 成n个子区间[xi,xi+1]，i=1,2,…,n，其中x1=a， xn+1=b。这样求定积分问题就分解为求和问题。
b ⎧ I = ∫ f ( x ' )dx ' ⎪ a ⎪ x ⎪ ' ' ⎨ I ( x) − I (a ) = ∫a f ( x )dx ⎪ ⎪ f ( x) = dI ( x) ⎪ dx ⎩
Xiamen University
Matlab and Engineering Calculation
Matlab and Engineering Calculation
2．高阶方法
MATLAB给出了quadl函数来替代quad8求定积分 [I,n]=quadl(fun,a,b,tol,trace) 其中参数的含义和quad函数相似。该函数可 以更精确地求出定积分的值，且一般情况下函数 调用的步数明显小于quad函数，从而保证能以更 高的效率求出所需的定积分值。
Xiamen University
Matlab and Engineering Calculation
例7-5 计算二重定积分 (1) 建立一个函数文件fxy.m： function f=fxy(x,y) global ki; ki=ki+1; %ki用于统计被积函数的调用次数 f=exp(-x.^2/2).*sin(x.^2+y); (2) 调用dblquad函数求解。 global ki;ki=0; I=dblquad('fxy',-2,2,-1,1) ki I= 1.57449318974494 ki = 1038
Trapezoidal numerical integration Z = trapz(X,Y) Z = cumtrapz(X,Y)
例7-4 已知函数在离散点的值，计算定积分值和积 分函数。 命令如下： x=pi/4:0.01:2*pi; y=sin(-x); I=trapz(x,y) Ix=cumtrapz(x,y) plot(x,y,’-o’,x,Ix)
矢量场旋度计算
[curlx,curly,curlz] = curl(X,Y,Z,U,V,W)
⎛ ∂Fz ∂Fy ⎞ ∇× F = i⎜ − ⎟+ ∂z ⎠ ⎝ ∂y ⎛ ∂Fx ∂Fz j⎜ − ∂x ⎝ ∂z ⎛ ∂Fy ∂Fx ⎞ ⎞ − ⎟ ⎟+k⎜ ∂x ∂y ⎠ ⎠ ⎝
Xiamen University