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 );
}