0

I was playing in Matlab when the following occured: I had an integral $I$ which I computed with the (a) Composite Simpson's Rule and (b) Adaptive Simpson's Method (my teacher told me that the built-in function quad in Matlab uses this method). These computations resulted in two approximations $I_1,I_2$ of $I$ for (a) and (b) respectively.

Prior to these computations I computed a reference solution to the problem, that is: an approximation of $I$ with a very low tolerance level (for example $10^{-15}$, compared to $10^{-6}$ which was used in later computations). Call this reference solution (approximation) $I_*$.

After finding $I_*$ I was able to compare $|I_1-I_*|$ with $|I_2-I_*|$. To my surprise $|I_1-I_*|<|I_2-I_*|$, and so: method (a) was more effective than method (b) in this particular example.

Now follows a general question: is there some general condition(s) for which it is true that $|I_1-I_*|<|I_2-I_*|$?

EDIT: Put

$y(x)=\frac{1-\arctan(p(x-1))/\pi}{2-\cos(\pi x)}$.

I am integrating $\pi \cdot y^2$ over $0 \leq x \leq 2.6$:

$$\pi \int_0^{2.6} y(x)^2 dx.$$

For $p=1$ I found that $|I_1-I_*|>|I_2-I_*|$, for $p=1000$ I found that $|I_1-I_*|<|I_2-I_*|$.

Below is the graph of $\pi \cdot y(x)^2$ for $p=1, 1000$

Graph of $\pi \cdot y(x)^2$ for $p=1, 1000$

bobbo
  • 59

1 Answers1

0

I suggest you have a look at

http://ezekiel.vancouver.wsu.edu/~cs330/lectures/integration/simpsons.pdf

They have a good presentation of what you ask.

Graeme
  • 3
Claude Leibovici
  • 260,315
  • This (see relation (33)) tells me (?) that given a tolerance level $\epsilon$, we should check that $|I_{[a,b]}-(I_{[a,c]}+I_{[c,b]})|<15\epsilon$ (where $I_{[a,b]}$ is the approximation of the integral over $[a,b]$) before proceeding. I do not understand how this answers my question, or hints at an answer to my question, or any other relation provided in the pdf-file. – bobbo Jan 08 '14 at 06:37
  • @bobbo. I agree : there is something strange I need to investigate. I shall come back soon. – Claude Leibovici Jan 08 '14 at 06:55
  • @bobbo. I do not use Mathlab but another CAS. I did not face any problem for the case where p=1. But, for p=1000, whatever I have been able to do, I have an error message claiming "Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small". So, I think that the problem is really related to your specific integral. Did you try the various options usable in "quad" ? Did you try using "quad1" and "quadgk" ? – Claude Leibovici Jan 08 '14 at 07:12
  • Yes, I also think it depends on the specific integrand. That is why there is such a change when we change $p$ from 1 to 1000. I can compute the integral with quad in Matlab, and using my own implementation of Composite Simpson's rule. What options do you mean? Here is another example (p. 5 of 6) for which $|I_1-I_|<|I_2-I_|$: http://www.math.usm.edu/lambers/mat460/fall09/lecture30.pdf ("Although the composite Simpson method contains half the error [!] of the adaptive quadrature method, 176 more function evaluations are required."). – bobbo Jan 08 '14 at 07:29
  • Also: there are som general remarks at p. 6 of 6 at http://www.math.usm.edu/lambers/mat460/fall09/lecture30.pdf . I do not understand them well. Maybe these remarks relates to my problem in some way (?). – bobbo Jan 08 '14 at 07:38
  • @bobbo. I think that the first and second remarks on page 6 are probably the keys for explaining what is going on with your problem. As I suggested, try your problem with plain Mathlab, invoking "quad", "quad1" and "quadgk". PLease, post the results as well as yours with a sufficient number of significant digits. I shall wait for that. Cheers. – Claude Leibovici Jan 08 '14 at 07:48
  • Did you plot Pi y^2 as a function of x for various values of p ? It looks terrible around x=1 for large values of p. May be you could split your [0,2.6] interval to [0,1]+[1,2.6]. Doing this, I do not face anymore my previous problem. – Claude Leibovici Jan 08 '14 at 08:15
  • Thanks! I will try that. I noticed that the relative effectiveness between Composite Simpson's Rule and Adaptive Simpson's Method decreases as $p$ goes from 1 to 1000. There was nothing wrong with my former code, but I used Composite Simpson's Method to find a reference solution. Now I computed a reference solution with quadl and I get these results: for $p=1$ and tolerance level 1e-6: abs(ValueQuad-reference),abs(ValueSim-reference)=0.000000381070319,0.003550437518000; for p=1000 (same tol.) I get 0.000001357413259, 0.001991390901978. – bobbo Jan 08 '14 at 08:26
  • Computing abs(ValueSimpson-reference)/abs(ValueQuad-reference) two times gives (for $p=1,1000$)

    1.0e+003 *

    9.317014047499411, 1.467048365144058.

    Thus, the relative effectiveness between Composite Simpson's Method and quad decreases (empirically) as $p$ goes from 1 to 1000.

     – bobbo Jan 08 '14 at 08:26