9--matlab数值微积分

数学实验

数值微积分

主要内容数值微(差)分 数值积分

数值定积分 数值重积分

数值差分1、数值差分

n维向量x=(x1, x2, , xn)的差分定义为 n-1维向量 x = (x2-x1, x3-x2, , xn- xn-1)。调用格式：Dx=diff(x) 如果x是向量，返回向量x的差分 ；如果x是矩阵，则按各列作差分。

DX=diff(A,n,dim):计算矩阵A的n阶 差分，dim=1时(缺省状态),按列计算 差分；dim=2,按行计算差分。 例如：

V=vander(1:6); V=V(:,end:-1:1) DV=diff(V) DV=diff(V,1,2) DV=diff(V,2) DV=diff(V,2,2)

%生成范德蒙行列式的矩阵 %按列计算V的一阶差分 %按行计算V的一阶差分 %按列计算V的二阶差分 %按行计算V的二阶差分

求下列定积分

1

01

e

x2

dx

1 1/2 = erf (1)* pi 2

sin( x) dx 0 x sin int( x) = x

= sin int(1)

0

sin(t ) dt t

数值积分求解定积分的数值方法多种多样，如 简单的梯形法、辛普生(Simpson) 等都是 法 经常采用的方法。它们的基本思想都是将 整个积分区间[a,b]分成n个子区间[xi,xi+1], i=1,2,…,n,其中x1=a,xn+1=b。这样求定 积分问题就分解为求和问题。

基本思想(1)分割: a=x0< x1< x2< < xn 1< xn =b, xi=xi xi 1; (2)近似代替: 小曲边梯形的面积近似为f(xi) xi (xi 1<xi<xi);A = lim (3)求和: 曲边梯形的面积近似为 f (x i ) x i ;n

(4)取极限: 设 =max{ x1, x2, , xn}, 曲边梯形的面积为A = lim f (x i ) x i 0 i =1n

0 i =1

Matlab 计算定积分函数介绍 Matlab 近似计算定积分的相关函数 1、梯形法：trapz

trapz(x,y)x 为分割点(节点)组成的向量， y 为被积函数在节点上的函数值组成的向量。

b

a

y y f ( x )dx = b a 0 y1 yn 1 n n 2 2

x = [ x0 , x1 , , xn ]

y = [ f ( x0 ), f ( x1 ), , f ( xn )]

trapz 举例 例：用梯形法计算下面定积分 ( 取 n=100)

dx 0 1 x 21

解： a=0, b=1, n=100, yi = f (xi) = 1/( 1+xi2 )做法思 想 trapz 函数

dx b a 0 1 x 2 n ( y0 / 2 y1 y2 yn 1 yn / 2)1

>> x=0:1/100:1; >> y=1./(1+x.^2); >> trapz(x, y)

trapz(x,1./(1+x.^2))

抛物线法2、 抛物线法：quad

b

a

f ( x)dx

quad(f,a,b,tol)f = f(x) 为被积函数，[a,b] 为积分区间，tol 为计算精度注： 不用自己分割积分区间 可以指定计算精度，若不指定，缺省精度是 10-6 精度越高，函数运行的时间越长 此处的函数 f 是数值形式，应该使用数组运算，即 点运算：.*,./ ,.\ ,.^将自变量看成是向量

quad 举例 例：用 quad 计算定积分：

dx 0 1 x 21

解：

>> quad('1./(1+x.^2)',0,1)>> quad('1./(1+x.^2)',0,1,10e-10) >> quad('1./(1+x.^2)',0,1,10e-16) 函数表达式一定要用 单引号 括起来！ 涉

及的运算一定要用 数组运算！

建立函数法求定积分例 求定积分：

3

0

e

0.5 x

sin( x )dx 6

(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)n为迭代次数。

洛巴托数值积分高精度Lobatto数值积分：quadl

调用格式： Q = quadl (FUN,a,b,TOL)不是数字1 其中参数的含义和quad函数相似

应用举例例：用函数trapz、quad和quadl分别数值积分 将计算结果精确值进行比较.

I=

1 2

/2

0

e

x2 2

dx

变限积分求导举例例：求 F '( x)

F ( x) = e sin(2 t )dtt 3 x

x2

syms x t F1=int(exp(t)*sin(2+sqrt(t^3)),x,0); F2=int(exp(t)*sin(2+sqrt(t^3)),0,x^2); Fi= F1+ F2; dF=diff(Fi)

综合应用举例例：求下列函数所围成的平面区域D的面积S, 并作出图形 ：

1 3 y = sin x, y = cos x, x = , x = 2 2x=-1/2:0.001:3/2;F1= sin(x); F2=cos(x); plot(x ,F1,'b-',x ,F2,'g-'), axis([-1,pi/4+1,-1.3,1.3]) xlabel('x'), ylabel('y'), fill([x;x],[F1;F2],'r') title(‘y=sinx , y=cosx和x=-0.5及x=1.5所围成的图形) syms x f1= cos(x)-sin(x); f2=-f1; S1 =int(f1,x,-0.5,pi/4); S2=int(f2,x, pi/4,1.5); S=S1+S2,Sj= double (S)