2016春计算方法实验指导-作业

第6次 数值积分

复化梯形公式：算法流程图

开始 定义f(x) 输入a, b, N (b-a)/N?h, 0 ?T 对k=1,2,…, N-1 T+f (a+k*h) ? T h* ?f (a)+2T+ f(b)? / 2 ? T 结束 复合辛普生公式：算法流程图

开始 定义f(x ) 输入a, b, N (b - a) / N?h, a+ h /2 ? x f(x) ? S1 , 0 ? S2 对k=1,2,…, N-1 S1 +f (a+k*h+h/2) ? S1 S2 +f (a+k*h) ? S2 h*? f (a)+4 S1 +2 S2+ f(b)? / 6 ? S2

结束

41

变步长梯形公式：算法流程图

开始 输入a,b,ε

b-a?h, h2?f(a)+f(b)??T1 0?S, a+h/2?x S+f(x)?S, x+h?x x

例6.1 给定积分

解：先输入主要初始参数

>>a=0.5; >>b=1;

>>f=inline('x^(1/2)');

%梯形公式

>>I1=(b-a)/2*(feval(f,a)+feval(f,b)) I1 =

0.426776695296637

%simpson公式

42

?10.5xdx,分别用梯形公式、Simpson公式、Cote公式作近似计算．

>>I2=(b-a)/6*(feval(f,a)+4*feval(f,(a+b)/2)+feval(f,b)) I2 =

0.430934033027025 %Cotes公式(n=4) >>tc=0;

>>C0=[7 32 12 32 7]; >>for i=0:4

tc=tc+C0(i+1)*feval(f,a+i*(b-a)/4); end

>>I3=(b-a)/90*tc I3 =

0.430964070495876 %准确值

>>I=int(char(f),a,b) >>vpa(I) I =

-1/6*2^(1/2)+2/3 ans =

0.43096440627115082519971854596505

例6.2 对积分I?sinx?4?0xdx,为使其精度达到10．

1若用复化梯形公式，应将[0,1]多少等分？ 若用复化Simpson公式，应将[0,1]多少等分？

解：直接按余项计算即可． 复化梯形公式的余项为:：

En(f)??复化Simpson公式余项为：

4b?a2''hf(?),12??(a,b)

b?a?h?(4)En(f)????(a,b),f(x)?C4[a,b] ??f(?),180?2?sinx对于f(x)?，在数学分析中，我们已证得以下不等式成立：

x1f(k)(x)?．

k?1直接利用上述不等式关系解答本题． 先输入误差精度：

43

>>Eps=1E-4 Eps =

1.000000000000000e-004

(1) 复化梯形公式

>>h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1))) % 先求出步长

h1 =

0.060000000000000

>>N1=ceil(1/h1) % 向上取整，得到等分区间数 N1 = 17

故可将区间17等分即可达到所要求的精度． (2) 复化Simpson公式

>>h2=power(Eps/abs(-(1-0)/180*1/(1+4)),1/4) %先求出步长 h2 =

0.547722557505166

>>N2=ceil(1/h2) % 向上取整，得到等分区间数 N2 = 2

故可将区间2等分即可达到所要求的精度．

? 扩展：

1) Matlab中复化梯形公式命令为I=trapz(x,y),复化Simpson公式命令为quad(). 2) Matlab中有四个取整函数,分别为ceil(),floor(),fix(),round(),分别表示向正无穷大方

向取整、向负无穷大方向取整、向靠近零方向舍入和四舍五入．

例6.3 对积分I?sinx?6dx??10,利用变步长方法求其近似值,使其精度达到． ?0x1解：利用变步长法前，先建立三种变步长复化积分公式的函数．注意在Matlab中直接用sin(0)/0得不到1，

sin(0)0matlab?NAN，因此解此题时，我们改用求极限的方法得到

函数值，此函数名为limit()．

先建立三种复化公式的函数文件，它们分别为复化梯形公式trap.m、复化Simpson公式为simpson.m、Cotes公式为cotes.m,三个函数的源程序如下： (1) 复化梯形公式trap.m

function T=trap(f,a,b,n) % 复化梯形公式求积分值 % f为积分函数 % [a,b]为积分区间

44