数值计算方法 实验报告4

学院：数学与计算机科学学院              专业：数学与应用数学

学号：                     姓名：

一．实验目的

1利用复化求积公式计算定积分，并比较误差；

2比较一阶导数和二阶导数的数值方法，并绘图观察特点.

二．实验题目

用复化梯形公式、复化辛普森公式、龙贝格公式求下列定积分，要求绝对误差为，并将计算结果与精度解进行比较：

   .

利用等距节点的函数值和端点的导数值，用不同的方法求下列函数的一阶和二阶导数，分析各种方法的有效性，并用绘图软件绘出函数的图形，观察其特点.

，  ， 

三．实验原理

1复化梯形公式

将积分区间剖分为等分，分点为，其中.在每个区间上用梯形公式，则有

            .

记.

2复化辛普森公式

将积分区间剖分为等分，分点为，其中.记区间的中点为，在每个区间上用辛普森公式，则得到所谓的复化辛普森公式：

，

即

.

3龙贝格公式的算法步骤为：

.输入及精度；

.置；

. 置，对分区间，并计算：

，；

.若不满足终止条件，做循环：

         ，

计算，

对计算：.

4向前差商公式：；

向后差商公式：；

中心差商公式：；

二阶导数公式：.

四．实验内容

实验一

第一小题：对于方程，利用程序shiyan1_01.m

内容如下：

%第一个函数的实验

clear

clc

fun=inline('(2/3)*x.^3.*exp(x.^2)');

S1=matrap(fun,1,2,170000);

S2=masimp(fun,1,2,250);

S3=maromb(fun,1,2,.5e-8);

s=exp(4);

Er1=abs(S1-s)

Er2=abs(S2-s)

Er3=abs(S3-s)

第二小题：对于方程，利用程序shiyan1_02.m

内容如下：

%第二个函数的实验

clear

clc

fun=inline('2*x./(x.^2-3)');

S1=matrap(fun,2,3,15000);

S2=masimp(fun,2,3,100);

S3=maromb(fun,2,3,.5e-8);

s=log(6);

Er1=abs(S1-s)

Er2=abs(S2-s)

Er3=abs(S3-s)

实验二

第一小题：对于方程，利用程序shiyan2_01.m

内容如下：

clear

clc

fun=inline('x.^5/20-(11./6)*x.^3');

dfun=inline('x.^4/4-(11./2)*x.^2');

ddfun=inline('x.^3-11*x');

n=8;h=2/n;

x=0:h:2;x1=x(2:n);

y=feval(fun,x);

dy=feval(dfun,x1);

ddy=feval(ddfun,x1);

for i=2:n

    dy1(i)=(y(i+1)-y(i))/h;

    dy2(i)=(y(i)-y(i-1))/h;

    dy3(i)=(y(i+1)-y(i-1))/(2*h);

    ddy1(i)=(y(i+1)-2*y(i)+y(i-1))/(h*h);

end

for i=1:n-1

    err1(i)=abs(dy1(i)-dy(i));

    err2(i)=abs(dy2(i)-dy(i));

    err3(i)=abs(dy3(i)-dy(i));

    errd2(i)=abs(ddy1(i)-ddy(i));  

end

[err1' err2' err3' errd2']

plot(x,y,'r')

hold on

plot(x1,dy,'y')

plot(x1,ddy,'k')

第二小题：对于方程，利用程序shiyan2_02.m

内容如下：

clear

clc

fun=inline('exp(-1./x)');

dfun=inline('(-1./x).*exp(-1./x)');

ddfun=inline('(-1./(x.^2)).*exp(-1./x)+1./(x.^2)');

n=8;h=2/n;

x=-2.5:h:-0.5;x1=x(2:n);

y=feval(fun,x);

dy=feval(dfun,x1);

ddy=feval(ddfun,x1);

for i=2:n

    dy1(i)=(y(i+1)-y(i))/h;

    dy2(i)=(y(i)-y(i-1))/h;

    dy3(i)=(y(i+1)-y(i-1))/(2*h);

    ddy1(i)=(y(i+1)-2*y(i)+y(i-1))/(h*h);

end

for i=1:n-1

    err1(i)=abs(dy1(i)-dy(i));

    err2(i)=abs(dy2(i)-dy(i));

    err3(i)=abs(dy3(i)-dy(i));

    errd2(i)=abs(ddy1(i)-ddy(i));  

end

[err1' err2' err3' errd2']

plot(x,y,'r')

hold on

plot(x1,dy,'y')

plot(x1,ddy,'')

五．实验结果

实验一

第一小题

T =

  146.5012        0         0          0          0          0          0           0

   83.9243   63.0653       0          0          0          0          0           0

   62.6132   55.5095   55.0058       0          0          0          0           0

   56.6535   54.6669   54.6108   54.6045       0          0          0           0

   55.1154   54.6027   54.5984   54.5982   54.5982       0          0           0

   54.7277   54.5984   54.5982   54.5982   54.5982   54.5982       0           0

   54.6305   54.5982   54.5982   54.5982   54.5982   54.5982   54.5982        0

   54.6062   54.5982   54.5982   54.5982   54.5982   54.5982   54.5982   54.5982

Er1 =

  4.5922e-009

Er2 =

  4.8409e-009

Er3 =

  1.4211e-014

第二小题

T =

    2.5000        0         0          0          0          0          0           0

    2.0192    1.8590       0          0          0          0          0           0

    1.85    1.8022    1.7984       0          0          0          0           0

    1.8088    1.7929    1.7922    1.7921       0          0          0           0

    1.7961    1.7918    1.7918    1.7918    1.7918       0          0           0

    1.7928    1.7918    1.7918    1.7918    1.7918    1.7918       0           0

    1.7920    1.7918    1.7918    1.7918    1.7918    1.7918    1.7918        0

    1.7918    1.7918    1.7918    1.7918    1.7918    1.7918    1.7918    1.7918

Er1 =

  4.9383e-009

Er2 =

  4.0302e-009

Er3 =

  1.0132e-012

实验二

第一小题

ans =

    0.2196    0.2196    0.2196    2.1920

    0.3627    0.8003    0.5815    2.1480

    0.5711    1.4367    1.0039    2.0560

    0.7667    2.0411    1.4039    1.9160

    0.9447    2.5991    1.7719    1.7280

    1.1003    3.0963    2.0983    1.4920

    1.2287    3.5183    2.3735    1.2080

    1.3251    3.8507    2.5879    0.8760

1.3847    4.0791    2.7319    0.4960

第二小题

ans =

    0.6932    0.6932    0.6932    0.1105

    0.4680    0.5532    0.5106    0.5030

    0.5236    0.6555    0.55    0.7793

    0.5907    0.8102    0.7005    1.2991

    0.6692    1.0727    0.8709    2.3982

    0.7473    1.6071    1.1772    5.1572

0.7567    3.0873    1.9220   14.2888

六．实验结果分析

1.利用复化辛普森公式比利用复化梯形公式，所取的更小，当达到相同精度时，利用辛普森公式等分次数更小，减少计算次数.

2.若利用同一公式，所取的大小与题设给出的精度之间的关系：当越大时，与精度之间的误差越小；反之，当越小时，与精度之间的误差越大。

3.龙贝格公式，是通过适当的线性组合，把复化梯形公式的近似值组合成更高精度的积分近似值.

七．附件

1 复化梯形公式程序：

%matrap.m

function s=matrap(fun,a,b,n)

%通途：用复化梯形公式求积分。

%格式：s=matrap(fun,a,b,n)  fun为被积函数，a,b为积分区间的左右端点，

% n为区间的等分数，s返回数值积分值

h=(b-a)/n;

s=0;

for k=1:n-1

    x=a+h*k;

    s=s+feval(fun,x);

end

s=h*(feval(fun,a)+feval(fun,b))/2.0+h*s;

2 复化辛普森公式程序：

%masimp.m

function s=masimp(fun,a,b,n)

%通途：用复辛普生形公式求积分。

%格式：s=masimp(fun,a,b,n)  fun为被积函数，a,b为积分区间的左右端点，

% n为区间的等分数，s返回数值积分值

h=(b-a)/n;

s1=0; s2=0;

for k=1:(n-1)

    x=a+h*k;

    s1=s1+feval(fun,x);

end

for k=0:(n-1)

    x=a+h*(k+1/2);

    s2=s2+feval(fun,x);

end

s=h/6*(feval(fun,a)+feval(fun,b)+2*s1+4*s2);

3 龙贝格公式程序：

%maromb.m

function s=maromb(fun,a,b,tol)

%用途：用龙贝格公式求积分

%格式：s=maromb(fun,a,b,tol), fun是被积函数，a, b是积分下、上限

%tol允许误差，s返回积分近似值

if nargin<4,tol=1e-4;end

i=1;j=1;h=b-a;

T(1,1)=h*(feval(fun,a)+feval(fun,b))/2;

T(i+1,1)=T(i,1)/2+sum(feval(fun,a+h/2:h:b-h/2))*h/2;

T(i+1,j+1)=(4^j*T(i+1,j)-T(i,j))/(4^j-1);

while abs(T(i+1,i+1)-T(i,i))>tol

    i=i+1;h=h/2; 

    T(i+1,1)=T(i,1)/2+sum(feval(fun,a+h/2:h:b-h/2))*h/2;

    for j=1:i

        T(i+1,j+1)=(4^j*T(i+1,j)-T(i,j))/(4^j-1);

    end    

end

T

s=T(i+1,j+1);