# Ninth Week Examples --- Quadrature:   Riemann Sums: Left Endpoint Rule; Midpoint Rule.  Trapezoid Rule. Simpson's Rule.

Following D. Milicic Tenth Week Examples.

# Numerical Integration

## Riemann Sums: Left endpoint, Right endpoint and Midpoint Rules.

```/* Milicic

Riemann sum--Left endpoint Rule

riemann1.c                           */

#include <stdio.h>

double f( double x) {

return(x*x) ;
}

int main() {

double a, b, delta, integral;
int i,n;

a = 0;
b = 1;
n = 100000;

delta = (b-a)/n;
integral = 0.0;

for (i=0; i < n; i=i+1)
integral = integral + f(a + i*delta)*delta;

printf("The integral of f from %f to %f is equal to %f\n",a,b,integral);
}

```

```
/* Milicic

Riemann sum--Right endpoint rule

riemann2.c                             */

#include <stdio.h>

double f( double x) {

return(x*x) ;
}

int main() {

double a, b, delta, integral;
int i,n;

a = 0;
b = 1;
n = 100000;

delta = (b-a)/n;
integral = 0.0;

for (i=0; i < n; i=i+1)
integral = integral + f(a + (i+1)*delta)*delta;

printf("The integral of f from %f to %f is equal to %f\n",a,b,integral);
}

```

## Midpoint Rule

/***************************************************************************** Treibergs Feb. 27, 2006 Midpoint Rule. Input endpoints a, b. Add up the areas of the strips with heights taken at the midpoints of the intervals. midpoint.c *****************************************************************************/ # include <stdio.h> # include <stdlib.h> # include <math.h> double f ( double x ); int main ( void ) { double a, am2, b, c, h, s; int i=1, n; printf ( "Midpoint Rule\n\n" ); printf ( " Enter the left and right endpoints : " ); scanf ( "%lf %lf", &a, &b ); if( a > b) /* swap a and b if a > b */ { c=a; a=b; b=c; } printf("\n Number of intervals Midpoint Approximation of Integral\n"); for ( n=1; n <= 20; n++) { h = (b - a) / n; am2 = a + h/2.0; s=0.0; for( i = 0; i < n; i++ ) s = s + f ( am2 + i*h ); printf (" %14d\t\t %21.15f\n", i, s*h); } printf ( "\t Actual int =%22.15f\n", 4.0*( atan(b) - atan(a) ) ); return EXIT_SUCCESS; } /* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */ double f ( double x ) { return 4.0/(1.0 + x*x) ; }

# Trapezoid Rule

```/*****************************************************************************
Treibergs                                                        Feb. 27, 2006

Trapezoid Rule. Input endpoints a, b. Integrate given function using trapezoid
rule with  n  intervals.

trapezoid1.c
*****************************************************************************/

# include <stdio.h>
# include <stdlib.h>
# include <math.h>

double
f ( double x );

int
main ( void )
{
double a, am2, b, c, h, h2, s;
int i=1, n;

printf ( "Trapezoid Rule\n\n" );
printf ( " Enter the left and right endpoints : " );
scanf ( "%lf %lf", &a, &b );

if( a > b)                         /*  swap a  and  b  if a > b   */
{
c=a;
a=b;
b=c;
}

printf("\n Number of intervals    Trapezoid Approximation of Integral\n");

for ( n=1;  n <= 20; n++)
{
h = (b - a) / n;
h2 = h / 2.0;
s=0.0;

for( i = 1; i < n; i++ )
s = s + f ( a + i*h );

printf (" %14d\t\t %21.15f\n", i, (f(a)+f(b))*h2+ s*h);
}
printf ( "\t    Actual int =%22.15f\n", 4.0*( atan(b) - atan(a) ) );

return EXIT_SUCCESS;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

double
f ( double x )
{
return 4.0/(1.0 + x*x) ;
}

```

```/*Milicic

trapezoid Rule

trapezoid.c          */

#include <stdio.h>

double f( double x) {

return(x*x) ;
}

int main() {

double a, b, delta, integral;
int i,n;

a = 0;
b = 1;
n = 100000;

delta = (b-a)/n;
integral = 0.0;

for (i=0; i < n; i=i+1)
integral = integral + (f(a + i*delta)+f(a+(i+1)*delta))/2*delta;

printf("The integral of f from %f to %f is equal to %f\n",a,b,integral);
}
```

# Simpson's Rule

```/*****************************************************************************
Treibergs                                                         Mar. 1, 2006

Simpson's Rule.

Input endpoints a, b. Integrate given function using Simpson's rule. Double
the number of intervals each step. Note that at each step, we only need
to evaluate the function at the points midway between the points from the
previous step. But this is nothing more than the midpoint rule. Then we use
the fact that Simpson's rule for 2n intervals is (trapezoid + 2 midpoint)/3.

simpson.c
*****************************************************************************/

# include <stdio.h>
# include <stdlib.h>
# include <math.h>

double
f ( double x );

int
main ( void )
{
double a, ap2, b, c, h, s, t;
int i = 1,  m, n = 1;

printf ( "Trapezoid / Midpoint / Simpson's Rules\n\n" );
printf ( " Enter the left and right endpoints : " );
scanf ( "%lf %lf", &a, &b );

if( a > b)                         /*  swap a  and  b  if a > b   */
{
c = a;
a = b;
b = c;
}

printf("\n       n      Trapezoid Approx.      Midpoint Approx.       Simpson's Rule\n");

h = b - a;
t = ( f ( a )  +  f ( b ) ) / 2.0;
for ( m = 1;  m <= 17; m++ )
{
ap2 = a  +  h / 2.0 ;
s = 0.0;

for( i = 0; i < n; i++ )
s = s  +  f ( ap2  +  i * h );

printf (" %9d %22.15f %22.15f %22.15f\n", n, t * h, s * h, (t  +  2.0 * s) * h / 3.0 );
t = t + s;
h = h / 2.0;
n = 2 * n;
}
printf ( "\t\t\t\t\t      Actual int =%21.15f\n", 4.0 * ( atan(b) - atan(a) ) );

return EXIT_SUCCESS;
}
/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

double
f( double x )
{
return 4.0 / ( 1.0  +  x * x );
}

/*  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =  =

Try the quadratic function like this instead! Simpson's rule is accurate on
quadratics by the way its defined. Be sure to replace the last printf by
printf ( "\t\t\t\t\t      Actual int =%21.15f\n", ( b * b * b  -  a * a * a ) / 3.0 );

double
f( double x )
{
return x*x;
}
*/
```

```/* Milicic

Simpson's method

simpson.c                */

#include <stdio.h>

double f(double x) {

return(x*x) ;
}

int main() {

double a, b, delta, integral;
int i,n;

a = 0;
b = 1;
n = 100000;

delta = (b-a)/n;
integral = 0.0;

for (i=0; i < n; i=i+1)
integral = integral + (f(a+i*delta)+4*f(a+(i+0.5)*delta)+f(a+(i+1)*delta))/6*delta;

printf("The integral of f from %lf to %lf is equal to %lf\n",a,b,integral);
}
```

```/*  Treibergs                  Mar. 6, 2006

Left-Right-Trapezoid Rules
These give lower, upper and trapezoid sum for increasing function.
We loop through the number of intervals for each sum. Then compute
the sum for that number of intervals. Note that the subtotal over the
middle  m - 1  intervals occurs in each of the three sums.

today.c                                   */

# include <stdio.h>
# include <stdlib.h>
# include <math.h>

double
f( double x );

int
main(void)
{
double s, a, b, d;
int  i, m, n = 20;

printf ( " Enter left and right endpoint : " );
scanf ( "%lf%lf", &a, &b );

printf ( " n \t\tLeft\t\t Right\t\t      Trapezoid\n" );

for( m = 2; m < = n; m++ )
{
s = 0.0;
d = ( b  -  a ) / m;
for ( i = 1; i < m;  i++ )
s = s + f ( a  +  i * d );

printf ( "%4d%21.15f%21.15f%21.15f\n", m, ( f ( a )  +  s ) * d, ( f ( b )  +  s ) * d,
( 0.5 * ( f ( a )  +  f ( b ) )  +  s ) * d  );
}

printf ( "  Actual integral = %47.15f\n", ( b * b * b  -  a * a * a ) / 3.0 );

return EXIT_SUCCESS;
}

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */

double
f( double x )
{
return  sqrt( .96 + x*x*x);
}

```

```/*****************************************************************************
Treibergs                                                         Mar. 8, 2006

Trapezoid and Midpoint Rule.

Input endpoints a, b. Integrate given function using trapezoid and midpoint
rules. Double the number of intervals each step. Note that at each step, we
only need to evaluate the function at the points midway between the points
from the previous step.

today1.c
*****************************************************************************/

# include <stdio.h>
# include <stdlib.h>
# include <math.h>

double
f ( double x );

int
main ( void )
{
double a, ahalf, b, c, e1, e2, h, s, t;
int i=1,  m, n=1;

printf ( " trapezRule\n\n" );
printf ( " Enter the left and right endpoints : " );
scanf ( "%lf %lf", &a, &b );

printf("\n  n \t Trapezoid Rule\t  Midpoint Rule\t       Midpt. Error    Trap. Error\n");

h = b  -  a;
e1 = h * h * h / 3.0;
e2 = 2.0 * e1;
t = ( f ( a )  + f ( b ) ) / 2.0;
for ( m = 1;  m <= 15; m++ )
{
ahalf = a  +  h / 2.0 ;
s = 0.0;

for( i = 0; i < n; i++ )
s = s + f ( ahalf  +  i * h );

printf ( "%5d%19.15f%19.15f  %14.6e  %14.6e\n", n, t*h, s*h, e1/n/n, e2/n/n);
t = t  +  s;
h = h / 2.0;
n = 2 * n;
}
printf ( "\t  Actual int =%21.15f\n", 4.0 * ( atan(b) - atan(a) ) );

return EXIT_SUCCESS;
}

double
f ( double x )
{
return 4.0 / ( 1.0  +  x * x );
}

```

```/*****************************************************************************
Treibergs                                                         Mar. 8, 2006
Simpson's Rule.

Input endpoints a, b. Integrate given function using Simpson's rule. Double
the number of intervals each step. Note that at each step, we only need
evaluate the function at the points midway between the points from the
previous step. Use the fact that Simpson's rule is (trap + 2 midpt)/3
For Simpson's rule with  n  intervals, the error is
E = (b-a)^5 M / ( 2880 n^4 )
where  M  = sup{ |f''''(c)| : c real }.  In case  f = 4/(1+x^2) then  M  = 96.
(Careful: in the calculus text "n" is something else: it is the number of
half intervals, which is our  2n.)

Recall that error estimates error of method and does not include round-off.

today.c
*****************************************************************************/
# include <stdio.h>
# include <stdlib.h>
# include <math.h>

double
f ( double x );

main ( void )
{
double a, ahalf, b, e, h, s, t;
i = 1, m, n = 1;

printf ( " Rule\n\n" );
printf ( " Enter the left and right endpoints : " );
scanf ( "%lf %lf", &a, &b );
printf ( "\n\tn \t  Simpson's Rule\t  Theoretical Error in Simpson's Rule\n" );

h = b  -  a;
e = h * h * h * h * h / 30.0;
t = ( f ( a )  +  f ( b ) ) / 2.0;

for ( m = 1;  m <= 15; m++ )
{
ahalf = a  +  h / 2.0 ;
s = 0.0;
for( i = 0; i < n; i++ )
s = s + f ( ahalf  +  i * h );

printf ("%11d %24.15f %28.15e\n", n, (t + 2*s) * h / 3, e / ( (double)n * n * n * n ) );

t = t  +  s;
h = h / 2.0;
n = 2 * n;
}

printf ( "Actual int =%24.15f\n", 4.0 * ( atan ( b ) - atan ( a ) ) );

return EXIT_SUCCESS;
}

/* * * * * * * * * * * * * * * * * * * * * * */

f( double x )
{
return 4.0 / ( 1.0  +  x * x );
}

```