Looks like you made a serious effort. First, let me find Boole, not something even most mathematicians would remember off hand. And I am one. :) Seriously, these higher order Newton-Cotes style rules are great for student assignments, and they were probably nice in the time of Boole. But you probably won't find a use for them ever again, except possibly in a homework assignment or on a test. What might you use? Actually, Gaussian quadrature rules are generally more useful, so you might want to look more closely at them when they come up in class. I've used various versions of Gaussian quadrature in many contexts over the years. Gaussian quadratires are significantly more efficient in terms of the number of function evals required for a given order rule.
Essentially, I see a weight out front of 2*h/45, with internal weights as the vector [7 32 12 32 7]. So a rule with a panel length of 5 points. As integration rules go, I suppose it is a decent one, with all positive weights, and none of them of a hugely different order of magnitude from the rest. (Those are good things in terms of noise amplification. However, if you were integrating noisy data, even this rule would amplify that noise, perhaps offsetting the reduction in integration error.)
Anyway, there is always an overlap at the joint of consecutive panels. So the composite rule would have weights like...
[7 32 12 32 14 32 12 32 14 32 12 32 14 32 12 32 14 ... 32 12 32 14 32 12 32 7]
Do you see the pattern? Simpson's rule works the same way. The first and last nodal weights end up being the normal value of 7, but the ones in between (at the join points) are doubled. That is because at the joint between panels, that same value is used twice.
For a composite rule with N panels here, you will have 4*N + 1 points in total. Thus one panel = 5 points. But two panels will use 9 points, 3 panels wil use 13 points, etc.
Now, how would I code it? Personally, I'd do it in a vectorized form. I'll show that first, doing the solution for a few different number of panels.
f = @(x) sin(x);
a = 0;
b = pi;
% remember, N is the number of panels
% save the result for various number of panels
int_est = zeros(5,1);
for N = 1:5
% compute the entire vector of weights at once.
% The first and last nodes get a different weight.
W = [7,repmat([32 12 32 14],1,N)];
W(end) = 7;
x = linspace(a,b,4*N + 1);
% the spacing between any pair of points
h = x(2) - x(1);
% the entire rule comes down to a dot product
int_est(N) = 2*h/45*dot(W,f(x));
end
N = (1:5)';
[N,4*N + 1,int_est]
ans =
1 5 1.99857073182384
2 9 1.99998313094599
3 13 1.99999858665236
4 17 1.99999975245457
5 21 1.99999993558359
As I recall, the true value for the integral over that interval is 2, so Boole does credibly well. Just to check how the gray haired memory is holding up these days...
syms X
int(sin(X),[0,pi])
ans =
2
Now, you want to use a loop? That too is not too difficult. Again, I'll do this for between 1 and 5 panels, but since you are wanting to use a loop over the panels, it will be a nested for loop. Note that there are hugely many ways you could write this, so I'll need to pick one.
f = @(x) sin(x);
a = 0;
b = pi;
% remember, N is the number of panels
int_est = zeros(5,1);
W = [7 32 12 32 7];
for N = 1:5
x = linspace(a,b,4*N + 1);
h = (b - a)/(4*N);
ind = (1:5);
% this inner loop goes over each panel,
% so things are still quasi-vectorized.
for m = 1:N
int_est(N) = int_est(N) + 2*h/45*sum(W.*f(x(ind)));
ind = ind + 4;
end
end
N = (1:5)';
[N,4*N + 1,int_est]
I won't display the result, since it is identical to the earlier more heavily vectorized version.
Note that I used a few different concepts here. In the first case, I used the function dot, in the second version, I used sum. I also computed h in two different ways, just for kicks.
Now, is the second version just as good as the first? Well, not perfectly so. Why not? because it evaluates the function TWICE at the join points between panels. So there are wasted function evals there. I suppose I could have been more intelligent and not done that. This is important SOME of the time, if a function eval is costly in time. When the few CPU cycles are irrelevant, I won't break a sweat.
Finally, where did you go wrong in what you did? Looking at the code you wrote, you were implicitly assuming that you need 4*N points in total. As you can see from what I wrote above, 4*N+1 is the correct number. At least that is the glaring error I saw in a brief glance.
xb = linspace(a,b,n/4);
==========================================================
Edit: I had the opportunity to see this answer, and I feel it was incomplete. How would i want to code a fully vectorized version of this rule? That is, If I were writing what I might consider to be a version of the code for my own use, what might I do to make it my own?
First, it needs to be a function. What I wrote were scripts. It needs to be well commented and fully documented. It needs to be vectorized. The best code would even test the function to see if it is vectorized, and if not, then would revert to evaluating the function in a loop. The best code is friendly code. It tries not to fail. If it does fail, it tries to inform the user about the problem, and perhaps even to suggest how to call it properly.
In a perfect world, I would write code that would be intelligent, choosing the number of panels in an intelligent way, or better yet, adaptively inserting new panels where the function was doing something interesting. In this case I might use a tiered version of the Booles rule, then using a Richardson extrapolation in a way that would allow me to both gain an implicitly higher order rule, as well as make the tool adaptive. The latter enhancements would be fun and interesting to write, but would potentially make for a very sophisticated code.
For here, I'll just make the code vectorized, so no loops will be used
function fint = myboole(f,a,b,N)
% myboole: Boole's rule numerical integration of a supplied function
% usage: fint = myboole(f,a,b,N)
%
% The total number of points that will be generated
% for evaluation by f will therefore be 4*N+1.
%
% arguments (input)
% f - a user supplied function handle.
%
% f is assumed to be evaluable in a vectorized form,
% where multiple points will be passed to it in one call.
% That requires the use of .* and ./ and .^ operators as
% necessary.
%
% a - scalar, numeric.
% The lower limit for the definite integral
%
% b - scalar, numeric.
% The upper limit for the definite integral
%
% N - (OPTIONAL) scalar, numeric, positive integer.
% The number of Boole's rule panels.
%
% Default: N == 1
%
% arguments (output)
% fint - the definite numerical integral over the
% interval [a,b], using N panels of Boole's rule.
%
% Boole's rule is a medium order Newton Cotes rule,
% It is closed, in the sense that the end points of
% the interval will be evaluated.
%
% https://en.wikipedia.org/wiki/Newton–Cotes_formulas
% https://en.wikipedia.org/wiki/Boole%27s_rule
%
% One panel of Boole's rule looks like:
% 2*h/45*(7*f(1) + 32*f(2) + 12*f(3) + 32*f(4) + 7*f(5))
%
% with an error term of
% -8/945*h^7*f''''''(zeta)
% where zeta is some point in the interval.
%
% Author: John D'Errico
% Date: 8/27/2020
% E-mail: woodchips@rochester.rr.com
% basic error checks
if nargin < 3
error('myboole:nargin','At least 3 input arguments MUST be provided')
end
if (nargin < 4) || isempty(N)
N = 1;
end
if ~isnumeric(a) || (numel(a) ~= 1)
error('myboole:invalida','a must be scalar, numeric')
end
if ~isnumeric(b) || (numel(b) ~= 1)
error('myboole:invalidb','b must be scalar, numeric')
end
% better code would be more flexible about the input function,
% for example, allowing a user supplied m-file. However
% an m-file can be trivially called as a function handle.
% better code will also verify if f is vectorized.
if ~isa(f,'function_handle')
error('myboole:invalidf','f must be a valid function handle')
end
% are the limits coincident? If so, then just return zero as the
% definite integral.
if a == b
fint = 0;
return
end
% N panels makes for 4*N+1 points in total
x = linspace(a,b,4*N+1);
% h will be negative if a > b. I could also have written
% h = (b - a)/(4*N);
h = x(2) - x(1);
% Nodal weights
W = [32 12 32 14];
W = [7,repmat(W,1,N)];
W(end) = 7;
% evaluate the function at the vector x
f_x = f(x);
% was f vectorized?
if size(f_x) ~= size(x)
if isvector(x) && (numel(f_x) == numel(x))
% f seems to have been vectorized, but swapped to a column vector
% be friendly and transpose it.
f_x = f_x.';
else
error('myboole:notvectorized','f must be a vectorized function')
end
end
% throw a warning message if nan or inf elements were found in f_x
if ~all(isfinite(f_x))
warning('myboole:finite','NaN or inf were encountered in the evaluation')
end
% form the dot product of the nodal weights and function values,
% scaling by 2*h/45. if a > b, then h will be negative, so will
% flip the sign on the result.
fint = sum(f_x.*W)*2*h/45;
end
Do you see anything interesting there? A huge feature of this code is it contains far more error checks, validity checks, etc., than it has effective code intended to compute a result. The comment lines are frequent, and the help is as long as I could make it, while still keeping it reasonable. in my eyes. Those comments will make the code more accessible to someone who wants to understand the code later, or if there is a bug found in the code one day.
Now, use it for varying values of the number of panels. I'll keep doubling the number of panels, which will effectively cut the step size h in half each time. Since we know the true value should be 2, I'll subtractt the result from 2 to plot it.
interr = zeros(1,10);
for ni = 1:10
interr(ni) = 2 - myboole(f,0,pi,2^(ni-1));
end
loglog(2.^(0:9),abs(interr))
xlabel '# of panels'
ylabel 'Absolute error in the integral estimate'
Do you see anything interesting there? As the number of panels doubles at each call, the step size h should be cut in half repeatedly. The integration error drops smoothlly and linearly on a loglog plot, up to a point.
However, after we reach 16 panels, that decrease tails off. It seems to completely stop beyond 32 panels. What happens here? This limit is imposed by double precision arithmetic.
I can use my own HPF toolbox to perform the computations in 100 significant digits.
DefaultNumberOfDigits 100
interr = hpf(zeros(1,10));
for ni = 1:10
interr(ni) = 2 - myboole(f,hpf(0),hpf('pi'),2^(ni-1));
end
loglog(2.^(0:9),abs(double(interr)),'r-o')
xlabel '# of panels'
ylabel 'Absolute error in the integral estimate'
Now performed in 100 decimal digits of precision, here you see the error decreases smoothly and linearly on the loglog plot. Eventually, that too would begin to fail with sufficient panels employed.
A little more analysis should allow us to investigate the error term, because the error should decrease with roughly the 7th power of the step size h.
Regardless, why did I add this final version? Because I want to teach people to write robust, friendly, efficient code. Comments are important, because they make your code readable.
