1

Using code I found from here for adaptive Simpson's method, Numerical Integration Via Adaptive Simpson's Method, I tried to implement it into my own equation, $(1 - (\sin^2[51/2])(\sin^2 x))^{-\frac12}$, but can't seem to get an output.

Simpson[f_, a_, b_, er_] := 
  iSimpson[f,(*N@*){a, (a + b)/2, b},(*N@*){f[a], f[(a + b)/2], f[b]},
    er];
iSimpson[f_, {a_, c_, b_}, {fa_, fc_, fb_}, er_] := 
  Module[{c1, c2, f1, f2, s1, s2, h2}, h2 = (b - a)/4; c1 = a + h2; 
   c2 = b - h2;
   f1 = f[c1];
   f2 = f[c2];
   s1 = 2 h2 (fa + 4 fc + fb)/3;
   s2 = h2 (fa + 4 f1 + 2 fc + 4 f2 + fb)/3;
   If[Abs[s1 - s2]/15 < er, s2 + (s2 - s1)/15, 
    iSimpson[f, {a, c1, c}, {fa, f1, fc}, er/2] + 
     iSimpson[f, {c, c2, b}, {fc, f2, fb}, er/2]]];

count = 0;
h[x_] := Module[{}, count++; (1 - (Sin[51/2]^2)*(Sin[x]^2))^(-0.5)];
N[Simpson[h, 0, Pi/2, 10.^-3] - 2]
count

it seems to work when I'm doing the method by hand, but can't get it work on my Mathematica.

Community
  • 1
johnson
  • 113
  • 4
  • Use Sin instead of sin; capitalization matters. – J. M.'s missing motivation Dec 15 '16 at 01:32
  • I changed both pi and sin to be capital and it's running a little bit but still can't get to the solution. – johnson Dec 15 '16 at 01:35
  • could it be because i'm plugging in 51 degrees at $Sin[51/2]$ but I'm doing an integral over radians? Nvm this wouldn't affect it – johnson Dec 15 '16 at 01:37
  • 1
    Read the documentation page for Sin to see that the argument should be in radians. You can use Sin[51/2 Degree] if you want to specify degrees. – Simon Rochester Dec 15 '16 at 01:41
  • I think it's working: Simpson[h, 0., Pi/2., 10.^-3] and NIntegrate[(1 - (Sin[51/2]^2)*(Sin[x]^2))^(-0.5), {x, 0, Pi/2}] agree to an error of less than 10^-3. (You subtract 2 from the integral in your second to last line, but I think that it is mistake. Is there a reason for it?) – Michael E2 Dec 15 '16 at 12:51

1 Answers1

0

Try it with a simple version!

f[x_] = (1 - Sin[51/2 Degree]^2 Sin[x]^2)^-0.5;
a = 0; b = π/2;
n = 10;
h = (b - a)/n;

simpson = 
 h/6 Sum[f[a + i h] + 4 f[a + (i + 1/2) h] + f[a + (i + 1) h], {i, 0, n - 1}] // N
{1.65231}

for comparison

NIntegrate[f[x], {x, 0, π/2}]
{1.65231}