数值积分与matlab求解

a
ab f ( x)dx (b a ) f ( ) 2
③ Simpson公式

b
a
1 ab f ( x)dx (b a) f (a) 4 f ( ) f (b) 6 2
在这三个公式中, 梯形公式把f(a), f(b)的加权平均值
1 f (a) f (b) 作为平均高度f()的近似值而获得的一种数 2
•
•
Fs= double (fi)
wQ1= double (abs(fi-Q1) ), wQ2= double (abs(fi-Q2) ) 值Q2、Q1和迭代次数FCNT14，取精度分别为104 ，Q2、Q1 分别与精 确值Fs 的绝对误差wQ2, wQ1
运行后屏幕显示I 的精确值Fs，用comsimpson.m 和quad.m 分别计算I 的近似
Sj =1.36203791318826
1.2
例2
变限积分的Matlab符号计算
x2 x
已知 F (x) e t sin(2 t 3 )dt ,求F ′(x)
解：输入程序： >> 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) 运行后屏幕显示计算变限积分F(x)的导数.
例3 用 MATLAB 的函数trapz 和 cumtrapz 分别计算


2 0
e - x sinxdx 精确到10-4 ，并与矩形公式比较.
解：将[0, π / 2]分成20 等份，步长为π / 40，输入程序如下（注意
trapz(y)是单位步长, trapz(y)*h=trapz(x,y)）： • >> h=pi/40; x=0:h:pi/2; y=exp(-x).*sin(x); • z1=sum(y(1:20))*h, z2=sum(y(2:21))*h, z=(z1+z2)/2 • z3=trapz(y)*h, z3h=trapz(x,y), z3c=cumtrapz(y)*h, • 运行后屏幕显示用矩形公式（9.3），（9.4）计算结果z1、z2 和 二者的平均数z、函数trapz和cumtrapz 分别计算结果z3、z3c.
• 例5 用 comsimpson.m 和quad.m 分别计算定积分
1 I ， 2
百度文库e
0
1 -
x2 2
dx
取精度为10-4 ，并与精确值比较.
解：输入程序 • • • • >> [Q1,FCNT14] = quad(@fun,0,1,1.e-4,3), Q2 =comsimpson (@fun,0,1,10000) syms x fi=int(exp( (-x.^2)./2)./(sqrt(2*pi)),x,0, 1);
1.3.1
梯形公式的Matlab程序
(一)用函数 trapz 计算定积分 • 调用格式一：Z =trapz(Y) • 调用格式二：Z =trapz(X,Y) • 调用格式三：Z = trapz (X,Y,DIM) 或 trapz (Y,DIM) (二) 用函数 cumtrapz 计算定积分 • 调用格式一：Z =cumtrapz (Y) • 调用格式二：Z =cumtrapz (X,Y) • 调用格式三：Z = cumtrapz (X,Y,DIM) 或 cumtrapz (Y,DIM)
按照这种思想，可构造出一些求积分值的近似公式。例
如 f ( ) 分别取
f ( ) f (
ab ) 2
和
f ( )
f (a) f (b) 2
则分别得到中矩形公式和梯形公式。 ① 梯形公式
1 f ( x)dx (b a) f (a) f (b) 2

b
a
②
矩形公式
b
1.3.3
辛普森公式的Matlab程序
• 调用格式一：quad(‘fun’,a,b) • 调用格式二：quad(‘fun’,a,b,tol) • 调用格式三：[Q,FCNT] = quad (...) • 调用格式四：quad(‘fun’,a,b, tol,TRACE) • 调用格式五：quad(‘fun’,a,b, tol,TRACE,P1,P2, …) 复合辛普森（Simpson）数值积分的MATLAB 主程序
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) 运行后屏幕显示计算面积的值 S 及其近似值Sj 如下 S =2*2^(1/2)+sin(1/2)-cos(1/2)-sin(3/2)-cos(3/2)
值积分方法。
矩形公式把[a,b] 的中点处函数值
f (
ab ) 2
作为平均高度f()的近似值而获得的一种数值积分方法。
Simpson公式是以函数f(x)在a, b, (a+b)/2这三点的函数值 f(a), f(b),
f( ab ) 2
的加权平均值 作为平均高度f()的近
似值而获得的一种数值积分方法。
复杂函数进行积分，这就是数值积分的思想，用代数插
值多项式去代替被积函数发f(x)进行积分是本章讨论数 值积分的主要内容。
1.1
定积分的Matlab符号计算
例1 由 y =sin x， y =cos x , x=-1/2,x=3/2所围成的平 面区域D.求平面区域D 的面积S.
解 输入作函数图形的程序 坐标调整 >> x=-1:0.001: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'), title('y=sinx , y=cosx 和x=-0.5及x=1.5所围成的平面区域的图形') 运行后屏幕显示图形. 求平面区域D 的面积S.输入计算面积S 的程序 >> syms x
计算


2 0
e - x sinxdx ，并与精确值比较.
解：将[0,/ 2]分成20 等份，步长为 / 40，输入程序如下（注意sum 和 cumsum的用法） • >> h=pi/40; x=0:h:pi/2; y=exp(-x).*sin(x); • z1=sum(y(1:20))*h,z2=sum(y(2:21))*h, • z=cumsum(y); z11=z(20)*h, z12=(z(21)-z(1))*h, 运行后屏幕显示计算结果分别如下 • z1 = z2 = z11 = z12 =0.3873 0.4036 0.3873 0.4036 求定积分的精确值，输入程序 • >> syms x • F=int(exp(-x)*sin(x),x,0, pi/2) • Fs= double (F) ,wz1=abs( Fs-z1), wz2= abs( Fs-z2) 运行后屏幕显示定积分的精确值Fs 和用矩形公式计算结果的绝对误差 wz1、wz2 分别如下 • F = Fs =1/2*(-1+exp(pi)^(1/2))/exp(pi)^(1/2) 0.3961 • wz1 = wz2 =0.0088 0.0075
• • • • • • • • • • • function y=comsimpson(fun,a,b,n) z1=feval (fun,a)+ feval (fun,b);m=n/2; h=(b-a)/(2*m); x=a; z2=0; z3=0; x2=0; x3=0; for k=2:2:2*m x2=x+k*h; z2= z2+2*feval (fun,x2); end for k=3:2:2*m x3=x+k*h; z3= z3+4*feval (fun,x3); end y=(z1+z2+z3)*h/3;
1 2 3 9 2 2 F ( x) x 2 x 3 x 2 x 3 ln( 2 x x 2 2 x 2 3 ) 4 16 16 2
(3) 被积函数f(x)没有具体的解析表达式, 其函数 关系由表格或图形表示。 对于这些情况, 要计算积分的准确值都是十分困难 的。由此可见, 通过原函数来计算积分有它的局限性, 因 而研究一种新的积分方法来解决Newton-Leibniz公式所不 能或很难解决的积分问题, 这时需要用数值解法来建立 积分的近似计算方法。 将积分区间细分,在每一个小区间内用简单函数代替
MATLAB求解数值积分
MATLAB求解连续函数积分
1 引言
Matlab求解连续函数积分
我们知道,若函数f(x)在区间[a,b]上连续且其原
函数为F(x),则可用Newton-Leibnitz公式

b
a
f ( x ) dx F (b) F ( a )
求定积分的值 , Newton-Leibnitz公式无论在理论上 还是在解决实际问题上都起了很大作用，但它并不 能完全解决定积分的计算问题，因为积分学涉及的 实际问题极为广泛，而且极其复杂，在实际计算中
1 ab ( f (a) 4 f ( ) f (b)) 6 2
（2）先用某个简单函数 ( x) 近似逼近f(x), 用 ( x) 代替原被积函数f(x)，即 b f ( x)dx b ( x)dx
a a
以此构造数值算法。从数值计算的角度考虑,函数 ( x) 应对f(x)有充分的逼近程度,并且容易计算其积分。由于多 项式能很好地逼近连续函数,且又容易计算积分,因此将 ( x) 选取为插值多项式, 这样f(x)的积分就可以用其插值多项式 的积分来近似代替
1.3
数值求积方法
建立数值积分公式的途径比较多, 其中最常用的有两 种：
（1）由积分中值定理可知，对于连续函数f(x)，在积分
区间[a,b]内存在一点ξ，使得 b f ( x)dx (b a) f ( )
a
a, b
即所求的曲边梯形的面积恰好等于底为(b-a)，高为 f ( ) 的矩形面积。但是点ξ的具体位置一般是未知的, 因 而 f ( ) 的值也是未知的, 称 f ( ) 为f(x) 在区间[a,b]上的平均 高度。那么只要对平均高度 f ( ) 提供一种算法，相应地就 获得一种数值求积方法
1.3.2
• • • • • •
矩形公式的Matlab程序
（一） 函数sum 的调用格式 调用格式一：sum(X) 调用格式二：sum (X,DIM) （二） 函数cumsum 的调用格式 调用格式一：cumsum(X) 调用格式二：cumsum (X,DIM)
例4 用 MATLAB 的函数sum 和 cumsum 及矩形公式
梯形数值积分的MATLAB 主程序 • function T=rctrap(fun,a,b,m) • n=1;h=b-a; T=zeros(1,m+1); x=a; • T(1)=h*(feval(fun,a)+feval(fun,b))/2; • for i=1:m • h=h/2; n=2*n; s=0; • for k=1:n/2 • x=a+h*(2*k-1); s=s+feval(fun,x); • end • T(i+1)=T(i)/2+h*s; • end • T=T(1:m);
MATLAB求解离散函数积分
2
Matlab求解离散函数积分
• • •
• •
1、S=trapz(X,Y) 等距梯形法求数值积分 调用格式一：Z =trapz(Y) 调用格式二：Z =trapz(X,Y) 调用格式三：Z = trapz (X,Y,DIM) 或 trapz (Y,DIM) 2、 S=cumsum(Y) 欧拉法求数值积分 调用格式一：cumsum(X) 调用格式二：cumsum (X,DIM) 3、 对于离散点的积分，先对其拟合获得函数表达式，再作为 连续被积函数求积分。
经常遇到以下三种情况：
(1) 被积函数f(x)并不一定能够找到用初等函数的 有限形式表示的原函数F(x)，例如：
1 sin x x2 dx和 e dx 0 x

1
0
Newton-Leibnitz公式就无能为力了 (2) 还有被积函数f(x)的原函数能用初等函数表示， 但表达式太复杂，例如函数 f ( x) x 2 2x 2 3 并不复杂，但积分后其表达式却很复杂，积分 后其原函数F(x)为：