常州大学数值分析第五章习题解答

化专012 闵建中 12102914

2012022302021342.(1)2()0

2()314==33

====A A A h h A A h A A h A A h A h ⎧⎪++=⎪--=⎨⎪⎪+=⎩

=解：

将f(x)=1,x,x 分别代入求积公式，并令左右相等，得

解得：，进一步，取f(x)x 代入上述求积公式，有左边右边0，但当f(x)x 时，上述积分公式左右不相等，因此该求积公式具有三次代数精度。 11221122.21

=1=*(13=2

311

==*(13)=

3311()((1)2(0)3())

3324

==,=1

39A A x x f x dx f f f -++-+≈-++⎰（） 解：

当f(x)时，令左右相等有：2）所以有当f(x)x 时，令左右相等有：0所以有所以求积公式为：当f(x)x 时，有左边右边，左右不相等，所以其代数精度为

4.解：

(1) : T=0.5*(3+1)=2 S=2 C=2 err1 = 0 err2=0

(2) : T=1*(1-3)=-2 S=-10/3 C=-10/3 err 1= 4/3 err2=0

(3) : T=0 S=0 C=0 err1=0 err2=0

(4) : T=e/2 S=1/6*(e+2exp(0.5)) C= 1 err1=-1.7183 err2= 0.0026

7.解：

复合梯形公式T 2n 的matlab 实现：

function I = trapezoid(fun,a,b,n)

n = 2*n;

h = (b-a)/(n-1);

x = a : h :b;

f = feval(fun , x);

I = h * (0.5*f(1)+sum(f(2:n-1))+ 0.5*f(n));

分别用复合梯形公式和复合辛普森公式计算积分，并与精确值比较，这里积分的

精确值由matlab int函数求出，不手工计算。

function trapezoid_and_sinpsom clc;

format long

syms x

Iexact = int(x*exp(x^2),x,0,1);

a=0;

b=1;

for n=2:1:4

t = trapezoid(@f,a,b,n)

s = simpson(@f,a,b,n)

err1 = vpa(Iexact - t,5)

err2 = vpa(Iexact - s,5)

end

function y=f(x)

y = x*exp(x^2);

return

从而得出第一小题的结果：

n = 2

t = 1.6697

s = 0.8695

err1 = -0.14144

err2 = -0.0018562

n = 3

t = 0.923798756293777

s = 0.859533825596209

err1 = -0.064658

err2 = -0.00039291

n = 4

t = 0.895892057505771

s = 0.859268455239111

err1 = -0.036751

err2 = -0.00012754

将f(x)简单修改就可以得到第2，3小题的结果8.解：

function test_romberg()

a=1;b=3;n=2;maxit=10;tol=0.5e-4; [R]=romberg(@f,a,b,n,tol,maxit);

R = R(1,4)

Rexact = log(3)

err = R-Rexact

return

function y=f(x)

y=1/x;

return

10.解：

高斯-勒让德法求解积分的matlab实现：

function I = guasslegendre(fun,a,b,n,tol)

% guasslegendre: guass-legendre method to calculation integral value

%Input : fun is the integrand input as a sring 'fun'

% a and b are upper and lower limits of integration % n is the initial number of subintervals

% tol is the tolerance

%Output : I is Integral result

% calculate quadrature node

syms x

p=sym2poly(diff((x^2-1)^(n+1),n+1))/(2^n*factorial(n));

tk=roots(p);

% calculate quadrature coefficient

Ak=zeros(n+1,1);

for i=1:n+1

xkt=tk;

xkt(i)=[];

pn=poly(xkt);

fp=@(x)polyval(pn,x)/polyval(pn,tk(i));

Ak(i)=quadl(fp,-1,1,tol);

end

% integral variable substitution, transformat [a, b] to [1, 1]

xk=(b-a)/2*tk+(b+a)/2;

% calculate the value of the integral function after variable substitution

fx=fun(xk)*(b-a)/2;

%calculation integral value

I=sum(Ak.*fx);

对第一小题，精确值由matlab int函数求出，不手工计算；

fun = inline('x.^5+3.*x.^2+2');