C++

Forum

 
Trapezoid method converges faster than Simpson method

KCl2000 (8)
Good afternoon,

I have been doing computer practices in C ++, and for an integration practice, the trapezoid method converges faster than the Simpson method. The function to be integrated is a first class elliptical integral of the form:

1/(\sqrt(1-k^2*sen^2theta))


Where k is bounded between [0,1), and the integration interval is between 0 and pi / 2. I have been thinking about the reason for this fact, and one of the hypotheses I have is that both the second derivative (which is used to determine the error level of the trapezoid method), and the fourth (which is used to determine the Simpson method error level), diverge in pi / 2, so it is not possible to calculate the maximum error level.

Does anyone know why this happens?
lastchance (6980)
KCl2000 wrote:
one of the hypotheses I have

What is the other one?


KCl2000 wrote:
both the second derivative ... and the fourth ... diverge in [at?] pi / 2

Since when? Not if k < 1. What value of k are you using?


Can we see your code? Are you running in single or double precision?


Also at
https://www.physicsforums.com/threads/trapezoid-method-converges-faster-than-the-simpson-method.982393/#post-6279251
Last edited on
KCl2000 (8)
Sorry, really meant, "This is the only hypothesis I have", my native language is not English and have some mistakes.

I am using several values of k (all of them less than 1), this phenomenon happens for all the values of k. The variables I am using are in all cases of type double.

This is the function for Simpson's method:

[
double intsimpson13(double a, double b, int n, double theta, double(*funcion)(double x, double y))
{

double h = (b-a)/double(n);

double sum1 = 0;
double sum2 = 0;
double x1,x2;
double fx1,fx2;

for (int i=1; i<(n/2); i++)
{
x1 = ((b-a)*2*i)/n+a;
fx1=funcion(x1,theta);
sum1=sum1+fx1;
}

for (int i=1; i<=(n/2); i++)
{
x2= ((b-a)*(2*i-1))/n+a;
fx2=funcion(x2,theta);
sum2=sum2+fx2;
}

double hint=2*sum1+4*sum2+funcion(a,theta)+funcion(b,theta);
return (hint*h)/3;
}][/code]

And this is the function I use to apply the trapezoid method:

[double inttrapecio(double a, double b, int n, double theta, double(*funcion)(double x, double y))
{

double h = (b-a)/double(n);

double sum = 0;
double x;
double fx;

for (int i=1; i<n; i++)
{
x = ((b-a)*i)/n;
fx=funcion(x,theta);
sum=sum+fx;
}

double inth = funcion(a,theta)+funcion(b,theta)+2*sum;
return inth*h/2;
}

][/code]

Thank you very much for the help,
Last edited on
lastchance (6980)
I believe that (give or take some bad algebra on my part) the second derivative is (read x for theta)

d2y/dx2=k2 cos(2x) (1-k2 sin2(x))-3/2 + (3/4)k4 sin2(2x)(1-k2 sin2(x))-5/2

and I could put a (very conservative) upper bound on that of
k2/(1-k2)5/2

So, as long as k is strictly less than 1 that won't diverge. The 4th derivative is similar.

I'll have a look at your code ... might take a while!
lastchance (6980)
1
2
 
x2= ((b-a)*(2*i-1))/n+a;
fx2=funcion(x2,theta);


Could you explain why you are calling a function of x and theta. Surely you should be calling a function of theta only? Or, at most, of (one of x or theta) and k. Are you muddling theta with k?

Can you show your code for the function that is to be integrated?
KCl2000 (8)
You are right, if k is less than 1, they do not diverge, but if it is observed that the fourth derivative has a much more drastic behavior than the second, I have represented them graphically but I don't know how to upload the file to the forum ...
KCl2000 (8)
Because k depends on the value of theta (it is sen (theta / 2)), the integration variable is x ...

The problem to solve is the integral of the period of a physical pendulum, and you have to modify theta from 0 to 140 degrees in intervals of five in five

Next I show the function to integrate:

[double funcion(double phi, double theta)
{
double k=sin(theta/2);
return 1/(sqrt(1-(k*k*sin(phi)*sin(phi))));
}][/code]
lastchance (6980)
KCl2000 wrote:
have represented them graphically but I don't know how to upload the file to the forum

OK, I am looking at the graph in your Physics Forum post at
https://www.physicsforums.com/threads/trapezoid-method-converges-faster-than-the-simpson-method.982393/#post-6279251

For the trapezium rule, error asymptotically proportional to n2, so
error = Cn-2
so
log(error)=constant-2 log(n)
so
slope is -2. yes, that is exactly the slope of your graph for the trapezium rule: so the trapezium rule is OK.

For Simpson's (1-3) rule, error asymptotically proportional to n4 (I think), so plotting log(error) against log(n) should give slope -4. Unfortunately, your graph does not have this slope.


As far as I can see, your code for Simpson's rule is correct - although it is unnecessarily complex; e.g.
x1 = ((b-a)*2*i)/n+a;
could be written as
x1 = 2.0*i*h+a;


The only things that I can see are:
(1) Some of your errors are actually VERY small (10-16 or so). It is possible that round-off error in more complicated calculations could exceed the theoretical truncation error on which the asymptotic behaviour is based.
(2) There are simply more arithmetic operations for Simpson's rule than the trapezium rule; maybe, again, the round-off is overwhelming the theoretical improvement with smaller intervals ("law of diminishing returns").


Basically, I can't see anything wrong with your code, other than that, from the slopes of the graphs, it seems to be Simpson's rule that is inconsistent rather than the trapezium rule. I can't see, obviously, what you have done to analyse your results, though I presume that you have checked this carefully.

EDIT - are you absolutely sure that you have the AXES THE RIGHT WAY ROUND in your Simpson's rule graphs. The Simpson 1-3 rule graph has a slope of -1/4 ... which is EXACTLY the reciprocal of -4 ... just saying!


One thing that you can do, though, is to apply Simpson's rule to a QUADRATIC function (parabola) - it should (theoretically) give an EXACT result, and any errors there would have to be associated with round-off and the finite accuracy with which floating-point numbers can be stored.



Last edited on
KCl2000 (8)
Thank you for your time and for the response,

In the first place, my main doubt is that, the formulas that limit the error in both the Simpson method and the trapeze method, establish the maximum value, but not the minimum value, then, it must necessarily converge in all Cases Simpson's method better than trapezoid?

On the other hand, how do you conclude that in Simpson's case the error is:

error = Cn-2

Wouldn't this be the maximum level, being able to take the error any value less than this?

In addition, the slope of the graph I think is approximately -0.5 in the case of the trapeze (instead of 2), and -0.25 in the case of Simpson's method.

On the other hand, as I have read, double type variables in C ++ can store between 15 and 16 decimal places, so that the error "stagnates" at approximately 10 ^ -16 is consistent, since arithmetic precision has been reached of the machine

I have checked Simpson's method for a quadratic function, and it gives me the exact method.

Thank you so much for everything
KCl2000 (8)
Unfortunately, I have put the axis well ...

That of these slopes is a coincidence for the value of k taken, if we take a value of different k other slopes come out (if k is small the function behaves better and both methods converge earlier).

If the forum let you upload files, you could upload the corresponding graph to other values of k ...
Last edited on
KCl2000 (8)

In the following links you can consult the graphics that I mention in my posts:

https://www.dropbox.com/s/uod51awznqq45g1/Graphics.docx?dl=0

https://www.dropbox.com/s/8h7nlzaoralxjdf/Graphics2.docx?dl=0
againtry (2313)
Imgur is a better place to display graphics than using docx. Chances are nobody will take it due to malware risk
lastchance (6980)
I think you are numerically integrating such a smooth, monotonic, well-behaved function that any general conclusion about trapezium rule versus Simpson's rule is fruitless. I've written a driver for your routines below. For the particular value of k chosen, round-off error overwhelms other errors by n=8 for trapezium rule.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
 
#include <iostream>
#include <iomanip>
#include <cmath>
using namespace std;

const double PI = 4.0 * atan( 1.0 );

// Routines supplied by KCL2000
double intsimpson13( double a, double b, int n, double theta, double(*funcion)(double x, double y));
double inttrapecio(double a, double b, int n, double theta, double(*funcion)(double x, double y));
double funcion(double phi, double theta);

// My own routines
double agm( double x, double y );
double ellipt( double k );
double mySimpson( double a, double b, int n, double theta, double (*f)( double, double ) );



int main()
{
   #define COL << setw( 25 ) << fixed << setprecision( 16 ) <<

   const double a = 0.0, b = PI / 2.0;
   const double theta = 70.0 * PI / 180.0;

   // Exact values
   double k = sin( 0.5 * theta );
   double exact = ellipt( k );
   cout COL "k="         COL k     << '\n';
   cout COL "ellipt(k)=" COL exact << "\n\n";

   // Values by numerical integration
   cout COL "n" COL "trapezium" COL "Simpson" COL "logE_trapezium" COL "logE_Simpson" << '\n';

   for ( int n = 4; n <= 60; n+= 2 )
   {
      double value_trapezium = inttrapecio ( a, b, n, theta, funcion );
      double value_Simpson   = intsimpson13( a, b, n, theta, funcion );
//    double value_Simpson   = mySimpson   ( a, b, n, theta, funcion );
      double logE_trapezium = log10( abs( value_trapezium - exact ) );
      double logE_Simpson   = log10( abs( value_Simpson   - exact ) );
      cout COL n COL value_trapezium COL value_Simpson COL logE_trapezium COL logE_Simpson << '\n';
   }
}


//****************** functions supplied by user KCL2000 **************

double intsimpson13( double a, double b, int n, double theta, double(*funcion)(double x, double y))
{
   double h = (b-a)/double(n);
   double sum1 = 0;
   double sum2 = 0;
   double x1,x2;
   double fx1,fx2;

   for (int i=1; i<(n/2); i++)
   {
      x1 = ((b-a)*2*i)/n+a;
      fx1=funcion(x1,theta);
      sum1=sum1+fx1;
   }

   for (int i=1; i<=(n/2); i++)
   {
      x2= ((b-a)*(2*i-1))/n+a;
      fx2=funcion(x2,theta);
      sum2=sum2+fx2;
   }

   double hint=2*sum1+4*sum2+funcion(a,theta)+funcion(b,theta);
   return (hint*h)/3;
}


double inttrapecio(double a, double b, int n, double theta, double(*funcion)(double x, double y))
{
   double h = (b-a)/double(n);
   double sum = 0;
   double x;
   double fx;

   for (int i=1; i<n; i++)
   {
      x = ((b-a)*i)/n;
      fx=funcion(x,theta);
      sum=sum+fx;
   }

   double inth = funcion(a,theta)+funcion(b,theta)+2*sum;
   return inth*h/2;
}


double funcion(double phi, double theta)
{
   double k=sin(theta/2);
   return 1/(sqrt(1-(k*k*sin(phi)*sin(phi))));
}

//****************** end of functions supplied by user KCL2000 **************



//****************** My routines ***********************

double agm( double x, double y )
{
   const double EPSILON = 1.0e-30;
   double a = x, g = y;
   double aold = a + 1, gold = g + 1;
   while ( abs( a - aold ) > EPSILON || abs( g - gold ) > EPSILON )
   {
      aold = a;
      gold = g;
      a = 0.5 * ( aold + gold );
      g = sqrt( aold * gold );
   }
   return a;
}


double ellipt( double k )
{
   return 0.5 * PI / agm( 1.0, sqrt( 1.0 - k * k ) );
}



double mySimpson( double a, double b, int n, double theta, double (*f)( double, double ) )
{
   double dx = ( b - a ) / n;
   double I = f( a, theta ) + f( b, theta );
   for ( int i = 1; i < n; i += 2 ) I += 4.0 * f( a + i * dx, theta );
   for ( int i = 2; i < n; i += 2 ) I += 2.0 * f( a + i * dx, theta );
   return I *= dx / 3.0;
}


                       k=       0.5735764363510460
               ellipt(k)=       1.7312451756570582

                        n                trapezium                  Simpson           logE_trapezium             logE_Simpson
                        4       1.7312451821583967       1.7312142958177228      -8.1869972160572022      -4.5103249679200763
                        6       1.7312451756575793       1.7312449237905148     -12.2830466849284292      -6.5988295180202945
                        8       1.7312451756570584       1.7312451734899457     -15.6535597745270216      -8.6641185449406333
                       10       1.7312451756570579       1.7312451756378409     -15.6535597745270216     -10.7163077561114104
                       12       1.7312451756570584       1.7312451756568843     -15.6535597745270216     -12.7597980124690782
                       14       1.7312451756570584       1.7312451756570566     -15.6535597745270216     -14.8084617345127647
                       16       1.7312451756570584       1.7312451756570582     -15.6535597745270216                     -inf
                       18       1.7312451756570579       1.7312451756570582     -15.6535597745270216                     -inf
                       20       1.7312451756570584       1.7312451756570582     -15.6535597745270216                     -inf
                       22       1.7312451756570588       1.7312451756570582     -15.1764385198073590                     -inf
                       24       1.7312451756570579       1.7312451756570582     -15.6535597745270216                     -inf
                       26       1.7312451756570582       1.7312451756570582                     -inf                     -inf
                       28       1.7312451756570584       1.7312451756570582     -15.6535597745270216                     -inf
                       30       1.7312451756570582       1.7312451756570582                     -inf                     -inf
                       32       1.7312451756570579       1.7312451756570582     -15.6535597745270216                     -inf
                       34       1.7312451756570582       1.7312451756570588                     -inf     -15.1764385198073590
                       36       1.7312451756570584       1.7312451756570579     -15.6535597745270216     -15.6535597745270216
                       38       1.7312451756570584       1.7312451756570584     -15.6535597745270216     -15.6535597745270216
                       40       1.7312451756570584       1.7312451756570582     -15.6535597745270216                     -inf
                       42       1.7312451756570586       1.7312451756570588     -15.3525297788630404     -15.1764385198073590
                       44       1.7312451756570588       1.7312451756570584     -15.1764385198073590     -15.6535597745270216
                       46       1.7312451756570588       1.7312451756570584     -15.1764385198073590     -15.6535597745270216
                       48       1.7312451756570582       1.7312451756570582                     -inf                     -inf
                       50       1.7312451756570582       1.7312451756570588                     -inf     -15.1764385198073590
                       52       1.7312451756570584       1.7312451756570584     -15.6535597745270216     -15.6535597745270216
                       54       1.7312451756570579       1.7312451756570584     -15.6535597745270216     -15.6535597745270216
                       56       1.7312451756570584       1.7312451756570582     -15.6535597745270216                     -inf
                       58       1.7312451756570579       1.7312451756570584     -15.6535597745270216     -15.6535597745270216
                       60       1.7312451756570582       1.7312451756570582                     -inf                     -inf

Last edited on
KCl2000 (8)
It gives you the same results as me, so I interpret that the results are correct ...

I have been reading, and for periodic functions, the trapeze method works very well due to cancellations that occur in the term of the error, I leave here some links to some files that I have found on the subject (there are some that did not get to understand them very well because I still don't master the math that is used):

https://www.math.drexel.edu/~tolya/301_periodic_integrands.pdf
https://math.mit.edu/~stevenj/trapezoidal.pdf
https://www.math.ucdavis.edu/~bremer/classes/fall2018/MAT128a/lecture7.pdf

Thank you very much for your help!!
againtry (2313)
Read 'Advanced Engineering Mathematics' 10th edition by Kreyszig - Wiley, Section 19.5 which discusses numerical integration error bounds, degree of precision and and numerical stability, and in doing so compares the rectangle, trapezoidal and Simpson rules.

Error bounds and there estimation gets back to the 'hypothesis' by OP - there is in fact a calculable upper and lower bound, and a separate calculation that can be made of roundoff errors and stability. Of course maybe a first class elliptical function doesn't fit the mold.
Topic archived. No new replies allowed.