实验10 数值积分
实验目的:
1.了解数值积分的基本原理; 2.熟练掌握数值积分的MATLAB 实现; 3.会用数值积分方法解决一些实际问题。
实验内容:
积分是数学中的一个基本概念,在实际问题中也有很广泛的应用。同微分一样,在《微积分》中,它也是通过极限定义的,由于实际问题中遇到的函数一般都以列表形式给出,所以常常不能用来直接进行积分。此外有些函数虽然有解析式,但其原函数不是初等函数,所以仍然得不到积分的精确值,如不定积分近似值,称为数值积分。
sin x
⎰ 0x d x 。这时我们一般考虑用数值方法计算其
1
10.1 数值微分简介
设函数y =f (x ) 在x 可导,则其导数为
*
f (x *+h ) -f (x *)
f '(x ) =lim (10.1)
h →0h
*
如果函数y =f (x ) 以列表形式给出(见表10-1),则其精确值无法求得,但可由下式求得其近似值
f (x *+h ) -f (x *)
f '(x ) ≈ (10.2)
h
*
表 10-1
一般的,步长h 越小,所得结果越精确。(10.2)式右端项的分子称为函数y =f (x ) 在
x *的差分,分母称为自变量在x *的差分,所以右端项又称为差商。数值微分即用差商近似
代替微商。常用的差商公式为:
f '(x 0) ≈
f (x 0+h ) -f (x 0-h )
(10.3)
2h
f '(x 0) ≈
-3y 0+4y 1-y 2
(10.4)
2h
f '(x n ) ≈
其误差均为O (h 2) ,称为统称三点公式。
y n -2-4y n -1+3y n
(10.5)
2h
10.2 数值微分的MATLAB 实现
MATLAB 提供了一个指令求解一阶向前差分,其使用格式为: dx=diff(x)
其中x 是n 维数组,dx 为n -1维数组[x 2-x 1, x 3-x 2, , x n -x 1], 这样基于两点的数值导数可通过指令diff(x)/h实现。对于三点公式,读者可参考例1的M 函数文件diff3.m 。
例1 用三点公式计算y =f (x ) 在x =1.0,1.2,1.4处的导数值,f (x ) 的值由下表给
解:建立三点公式的M 函数文件diff3.m 如下:
function f=diff3(x,y) n=length(x);h=x(2)-x(1); f(1)=(-3*y(1)+4*y(2)-y(3))/(2*h); for j=2:n-1
f(j)=(y(j+1)-y(j-1))/(2*h); end
f(n)=(y(n-2)-4*y(n-1)+3*y(n))/(2*h);
在MATLAB 指令窗中输入指令:
x=[1.0,1.1,1.2,1.3,1.4];y=[0.2500,0.2268,0.2066,0.1890,0.1736];diff3(x,y)
运行得各点的导数值为:-0.2470,-0.2170,-0.1890,-0.1650,-0.0014。所以y =f (x ) 在x =1.0,1.2,1.4处的导数值分别为-0.2470,-0.1890和-0.0014。
对于高阶导数,MATLAB 提供了几个指令借助于样条函数进行求导,详细使用步骤如下: step1:对给定数据点(x,y ),利用指令pp=spline(x,y),获得三次样条函数数据pp ,供后面ppval 等指令使用。其中,pp 是一个分段多项式所对应的行向量,它包含此多项式的阶数、段数、节点的横坐标值和各段多项式的系数。
step2:对于上面所求的数据向量pp ,利用指令[breaks,coefs,m,n]=unmkpp(pp)进行处理,