MATLAB: Can I make this code for Boole’s rule smaller

newton-cotes boole integration

I'm trying to write a script for the composite form of Boole's rule:
∫ g(t) dt = (2h/45)*[7g((a)+(a+4h))+32g((a+h)+(a+3h))+12g(a+2h)]
I have this code (below) that computes Boole's rule. I don't think I'm using the composite of Boole's rule though and I'd like to. Can anyone suggest how I could change my code? Any help would be much appreciated.
function r = boole_test(f,a,b,n)
h = (b - a) / (n * 4);
r = 7 * f(a);
x = a + h;
for i = 1 : n-1
    r = r + 32 * f(x);
    x = x + h;
    r = r + 12 * f(x);
    x = x + h;
    r = r + 32 * f(x);
    x = x + h;
    r = r + 14 * f(x);
    x = x + h;
end
r = r + 32 * f(x);
x = x + h;
r = r + 12 * f(x);
x = x + h;
r = r + 32 * f(x);
r = r + 7 * f(b);
r = r * h*2/45;

Best Answer

First of all, you need to actually know what Boole's Rule is!
http://mathworld.wolfram.com/BoolesRule.html
Why did I know I should check? Because there is a trivial test for a numerical integration rule. If it does not integrate a constant function correctly, you got something wrong!
ALWAYS do simple tests like this. ALWAYS.
The weights on your points are wrong. They must be, because if you put a unit constant function in there, thus f(x) == 1, over a paneled region [0,4*h], you should get
2*h/45*(90) = 4*h
So the weights must sum to 90. The weights that you had sum to 53. So even before I did anything, I knew you were wrong. As well, weights for a Newton-Cotes rule like this will always be symmetric. So another quick test failed.
Looking up the correct form of Boole's Rule gives us...
2*h/45*( 7*f(x) + 32*f(x+h) + 12*f(x+2*h) + 32*f(x+3*h) + 7*f(x+4*h) )
First test: Are they symmetric? Yes.
Second test:
sum([7 32 12 32 7])
ans =
    90
Basic tests. Common sense. Always apply common sense.
Finally, how do you form a composite Boole's rule?
Remember that at the junction of the panels, the 4*h nodes will always be used twice. This happens with Simpson's Rule, it even does with trapezoidal rule. So trapezoidal rule can be seen to have weights like this:
[1 2 2 2 2 ... 2 2 2 2 2 1]
Simpson's rule?
[ 1 4 2 4 2 ... 4 2 4 2 4 1]
The first and last points end up with half the weight that the internal points get. So even though a single panel of Simpsons rule has relative weights that look like weights of [1 4 1], in the composite rule, the joints between the panels get double weight.
The same applies to Boole's rule. Weights of [7 32 12 32 7] turn into a sequence like this:
[7 32 12 32 14 32 12 32 14 ,,, 32 12 32 14 32 12 32 7]
So the joints have a weight of 14. But the first and last points get a weight of 7. THINK ABOUT IT!
How do you implement it? Easy peasy. No explicit for loops required.
Iboole = 2*h/45*(14*sum(f(1:4:end)) + 32*sum(f(2:4:end)) + 12*sum(f(3:4:end)) + 32*sum(f(4:4:end)));
Iboole = Iboole - 2*h/45*(7*f(1) + 7*f(end));
That second line just corrects for the first and last points. Test it out. Integrate the function f(x) == 1, over the interval with limits [0,12], thus a stride of 1. The integral should be 12.
f = ones(1,13);
h = 1;
Iboole = 2*h/45*(14*sum(f(1:4:end)) + 32*sum(f(2:4:end)) + 12*sum(f(3:4:end)) + 32*sum(f(4:4:end)));
Iboole = Iboole - 2*h/45*(7*f(1) + 7*f(end));
Iboole
Iboole =
         12
Can you do this using a loop? Of course. I can even do it using the weights where I do not double the weights at the joints between panels. That is because here, I will use those joint nodes TWICE. Again, think about it. Look carefully at what was done here.
n = numel(f);
if mod(n,4) ~= 1
  error('The sky is falling! f is the wrong length')
end
Iboole = 0;
for i = 1:4:(n-4)
  Iboole = Iboole + 7*f(i) + 32*f(i+1) + 12*f(i+2) + 32*f(i+3) + 7*f(i+4);
end
Iboole = Iboole*2*h/45;
Again, when I tested it, I got:
Iboole
Iboole =
    12
Note that n MUST ALWAYS be 1 plus a multiple of 4. So 5, 9, 13, 17, etc.
Always apply common sense tests to your work. Use such tests at every step of the way. This helps to ensure that you don't get something crazy out at the end, and then need to track through everything to find the bug. Here the very first bug was you had the wrong coefficients. You would never have found that by simply trying to debug your code.
Related Question