数值计算课程PROJECT- 牛顿-柯特斯公式(Newton-Cotes Methods)

数值计算课程PROJECT

牛顿-柯特斯公式(Newton-Cotes Methods)

例题介绍：

请在MATLAB环境下编写牛顿-柯特斯公式，以求得π的近似值。

Instruction:

Using the Newton-Cotes method to solve all the following problems.

Solving the problems with MATLAB, printing out your MATLAB code, figures, as well as necessary problem-solving procedures on paper, and submit before final exam. (No late submission is accepted) Students must solve all these questions correctly to get 5 point extra credits, no partial credit is given. Plagiarizing from other’s work will be treated as 0.

Problems I: The value of π can be calculated from the integral 1

2121dx x π-=+⎰(a)Evaluate the integral by using rectangle method, using 60 subintervals.

(b)Evaluate the integral by using midpoint method, using 60 subintervals.

(c)Evaluate the integral by using trapezoidal method, using 60 subintervals.

(d)Evaluate the integral by using Simpson’s 1/3 method, using 60 subintervals.

(e)Evaluate the integral by using Simpson’s 3/8 method, using 60 subintervals.

(f)Compare the results and discuss the error from each method.

MATLAB CODE:

f = @(x) 2/(1 + x^2); n = 60; % #subintervals

a = -1; % min

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

xi = a:h:b;

part A: 长方形法则

la = 0;

for i=1:length(xi)-1

la = la + h*f(xi(i));

end

fprintf('A. Rectangle method approx: %f\n', la)

part B: 中点法则

lb = 0;

for i = 1:length(xi)-1

lb = lb + h*f((xi(i) + xi(i+1))/2);

end

fprintf('B. Midpoint method approx: %f\n', lb)

part C: 梯形法则

lc = 0;

for i=1:length(xi)-1

lc = lc + h*(f(xi(i)) + f(xi(i+1)))/2;

end

fprintf('C. Trap. method approx: %f\n', lc)

part D: 辛普森1/3法则

ld = (h/3)*(f(xi(1)) + f(xi(end)));

for i = 2:2:length(xi)-1

ld = ld + (h/3)*4*f(xi(i));

end

for i = 3:2:length(xi)-2

ld = ld + (h/3)*2*f(xi(i));

end

fprintf('D. Simpsons 1/3 method approx: %f\n', ld)

part E: 辛普森3/8法则

le = (3*h/8)*(f(xi(1)) + f(xi(end)));

for i = 2:3:length(xi)-1

le = le + (3*h/8)*3*(f(xi(i)) + f(xi(i+1)));

end

for i = 3:3:length(xi)-3

le = le + (3*h/8)*2*f(xi(i));

end

fprintf('E. Simpsons 3/8 approx: %f\n', le)

part F: 误差分析

fprintf('\nF.')

fprintf('Relative error for A is %f %%\n', abs(la-pi)*100/pi) fprintf('Relative error for B is %f %%\n', abs(lb-pi)*100/pi) fprintf('Relative error for C is %f %%\n', abs(lc-pi)*100/pi) fprintf('Relative error for D is %f %%\n', abs(ld-pi)*100/pi) fprintf('Relative error for E is %f %%\n', abs(le-pi)*100/pi)

--------------------------------------------------------

结果如下：(Command Window results below)

A. Rectangle method approx: 3.141407

B. Midpoint method approx: 3.141685

C. Trap. method approx: 3.141407

D. Simpsons 1/3 method approx： 3.141593

E. Simpsons 3/8 approx: 3.141301

F.

Relative error for Rectangle method is 0.005895

Relative error for Midpoint method is 0.002947

Relative error for Trapezoidal method is 0.005895

Relative error for Simpsons 1/3 is 0.000000

Relative error for Simpsons 3/8 is 0.009284