数值分析(第五版)计算实习题第四章作业

1、（1）：复合梯形

建立m文件：

function t=natrapz(fname,a,b,n) 

h=(b-a)/n;  

fa=feval(fname,a);fb=feval(fname,b);

f=feval(fname,a+h:h:b-h+0.001*h); 

t=h*(0.5*(fa+fb)+sum(f)); 

输入：

>> syms x

>> f=inline('sqrt(x).*log(x);');

>> natrapz(f,eps,1,10)

输出：

ans =

  -0.417062831779470

输入：

>> syms x

>> f=inline('sqrt(x).*log(x);');

>> natrapz(f,eps,1,100)

输出：

ans =

  -0.443117908008157

输入：

>> syms x

>> f=inline('sqrt(x).*log(x);');

>> natrapz(f,eps,1,1000)

输出：

ans =

  -0.4443875397162

复合辛普森

建立m文件：

function t=comsimpson(fname,a,b,n)

h=(b-a)/n;

fa=feval(fname,a);fb=feval(fname,b);

f1=feval(fname,a+h:h:b-h+0.001*h);

f2=feval(fname,a+h/2:h:b-h+0.001*h); 

t=h/6*(fa+fb+2*sum(f1)+4*sum(f2));

输入：

>> syms x

>> f=inline('sqrt(x).*log(x);');

>> format long；

>>comsimpson(f,eps,1,10)

输出：

ans =

  -0.43529700746

输入：

>>syms x

>>f=inline('sqrt(x).*log(x);');

>>comsimpson(f,eps,1,100)

输出：

ans =

  -0.444161178415673

输入：

>>syms x

>>f=inline('sqrt(x).*log(x);');

>>comsimpson(f,eps,1,1000)

输出：

ans =

  -0.444434117614180 

（2）龙贝格

建立m文件：

function [RT,R,wugu,h]=Romberg(fun,a,b,wucha,m) 

%RT是龙贝格积分表

%R是数值积分值

%wugu是误差估计

%h是最小步长

%fun是被积函数

%a b是积分下、上限

%m是龙贝格积分表中行最大数目

%wucha是两次相邻迭代值的绝对误差限

n=1;h=b-a;wugu=1;x=a;k=0;RT=zeros(4,4);

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

while((wugu>wucha)&(k    k=k+1;h=h/2;s=0;

    for j=1:n

        x=a+h*(2*j-1);s=s+feval(fun,x);

    end

    RT(k+1,1)=RT(k,1)/2+h*s;n=2*n;

    for i=1:k

        RT(k+1,i+1)=((4^i)*RT(k+1,i)-RT(k,i))/(4^i-1);

    end

    wugu=abs(RT(k+1,k)-RT(k+1,k+1));

end

R=RT(k+1,k+1);

输入：

>>fun=inline('sqrt(x).*log(x)');

>> [RT,R,wugu,h]=Romberg(fun,eps,1,1e-5,13)

输出：

RT =

  1 至 5 列

  -0.000000268546145                   0                   0                   0                   0

  -0.2450670140209  -0.3267528040047                   0                   0                   0

  -0.358104125949240  -0.395783944552250  -0.400386020588741                   0                   0

  -0.408090073087781  -0.424752055467295  -0.426683262861631  -0.4271006794055                   0

  -0.429474601629505  -0.436602777810080  -0.437392825966266  -0.437562819031419  -0.437603847029951

  -0.4383494461832  -0.441361125405941  -0.4416783485799  -0.441746372747455  -0.441762778840459

  6 列

                   0

                   0

                   0

                   0

                   0

  -0.441766844267449

R =

  -0.441766844267449

wugu =

     4.0654269774412e-06

h =

.0312********

（3）自适应辛普森

输入：

>> f=inline('sqrt(x).*log(x)');

>> q=quad(f,0,1,1e-4)

输出：

q =

  -0.443975572951728

2.（1）

复合辛普森

建立m文件

function q=combinesimpson2(F,x0,a,b,n)

%复合Simpson多元求积公式

%F—被积函数

%x0—被积函数自变量

%[a,b]积分区间

%n—区间份数

x=linspace(a,b,n+1);

q=0;

for k=1:n

    q=q+subs(F,x0,x(k))+4*subs(F,x0,(x(k)+x(k+1))/2)+subs(F,x0,x(k+1));

end

q=q*(b-a)/n/6;

输入：

>> clear

>> syms x y;

>> F=exp(-x.*y);

>> s=combinesimpson2(combinesimpson2(F,'x',0,1,4),'y',0,1,4)

输出：

s =

exp(-1)/576 + exp(-1/2)/144 + exp(-1/4)/72 + exp(-3/4)/144 + exp(-1/8)/36 + exp(-3/8)/36 + exp(-5/8)/72 + exp(-7/8)/72 + (5*exp(-1/16))/144 + exp(-3/16)/24 + exp(-5/16)/36 + exp(-7/16)/36 + exp(-9/16)/144 + exp(-1/32)/36 + exp(-3/32)/18 + exp(-5/32)/36 + exp(-7/32)/36 + exp(-9/32)/36 + exp(-15/32)/36 + exp(-21/32)/36 + exp(-1/)/36 + exp(-3/)/18 + exp(-5/)/18 + exp(-7/)/18 + exp(-9/)/36 + exp(-15/)/18 + exp(-21/)/18 + exp(-25/)/36 + exp(-35/)/18 + exp(-49/)/36 + 47/576

>> double(s)

ans =

   0.796599967946203

高斯求积公式

function q=gaussquad(F,x0,a,b,n)

%Gauss求积公式

%F—被积函数

%x0—被积函数自变量

%[a,b]积分区间

%n—节点个数

syms t;

F=subs(F,x0,(b-a)/2*t+(a+b)/2);

[x,A]=gausspoints(n);

q=(b-a)/2*sum(A.*subs(F,t,x));

输入：

>> clear

>> syms x y;F=exp(-x.*y);

>> s=gaussquad(gaussquad(F,x,0,1,4),y,0,1,4)

输出：

s =

0.7966

（2）复合辛普森

输入：

>> syms x y;

>> f=exp(-x.*y);

>> s=combinesimpson2(combinesimpson2(f,y,0,sqrt(1-x^2),4),x,0,1,4)

 输出：

s =

(3^(1/2)*(exp(-3^(1/2)/4) + 2*exp(-3^(1/2)/8) + 2*exp(-3^(1/2)/16) + 2*exp(-(3*3^(1/2))/16) + 4*exp(-3^(1/2)/32) + 4*exp(-(3*3^(1/2))/32) + 4*exp(-(5*3^(1/2))/32) + 4*exp(-(7*3^(1/2))/32) + 1))/576 + (7^(1/2)*(exp(-(3*7^(1/2))/16) + 2*exp(-(3*7^(1/2))/32) + 2*exp(-(3*7^(1/2))/) + 2*exp(-(9*7^(1/2))/) + 4*exp(-(3*7^(1/2))/128) + 4*exp(-(9*7^(1/2))/128) + 4*exp(-(15*7^(1/2))/128) + 4*exp(-(21*7^(1/2))/128) + 1))/1152 + (15^(1/2)*(exp(-15^(1/2)/16) + 2*exp(-15^(1/2)/32) + 2*exp(-15^(1/2)/) + 2*exp(-(3*15^(1/2))/) + 4*exp(-15^(1/2)/128) + 4*exp(-(3*15^(1/2))/128) + 4*exp(-(5*15^(1/2))/128) + 4*exp(-(7*15^(1/2))/128) + 1))/1152 + (15^(1/2)*(exp(-(7*15^(1/2))/) + 2*exp(-(7*15^(1/2))/128) + 2*exp(-(7*15^(1/2))/256) + 2*exp(-(21*15^(1/2))/256) + 4*exp(-(7*15^(1/2))/512) + 4*exp(-(21*15^(1/2))/512) + 4*exp(-(35*15^(1/2))/512) + 4*exp(-(49*15^(1/2))/512) + 1))/1152 + (39^(1/2)*(exp(-(5*39^(1/2))/) + 2*exp(-(5*39^(1/2))/128) + 2*exp(-(5*39^(1/2))/256) + 2*exp(-(15*39^(1/2))/256) + 4*exp(-(5*39^(1/2))/512) + 4*exp(-(15*39^(1/2))/512) + 4*exp(-(25*39^(1/2))/512) + 4*exp(-(35*39^(1/2))/512) + 1))/1152 + (55^(1/2)*(exp(-(3*55^(1/2))/) + 2*exp(-(3*55^(1/2))/128) + 2*exp(-(3*55^(1/2))/256) + 2*exp(-(9*55^(1/2))/256) + 4*exp(-(3*55^(1/2))/512) + 4*exp(-(9*55^(1/2))/512) + 4*exp(-(15*55^(1/2))/512) + 4*exp(-(21*55^(1/2))/512) + 1))/1152 + (63^(1/2)*(exp(-63^(1/2)/) + 2*exp(-63^(1/2)/128) + 2*exp(-63^(1/2)/256) + 2*exp(-(3*63^(1/2))/256) + 4*exp(-63^(1/2)/512) + 4*exp(-(3*63^(1/2))/512) + 4*exp(-(5*63^(1/2))/512) + 4*exp(-(7*63^(1/2))/512) + 1))/1152 + 1/24

>> double(s)

ans =