數值積分方法 長方形(左邊 , 右邊 , 中央) 梯形法 , 辛普森 方法

Write functions to calculate the definite integral of a function ƒ(x) using all five of the following methods:

• rectangular
  • left
  • right
  • midpoint
• trapezium
• Simpson's
  • composite

Your functions should take in the upper and lower bounds (a and b), and the number of approximations to make in that range (n).

Demonstrate your function by showing the results for:

• ƒ(x) = x3, where x is [0,1], with 100 approximations. The exact result is 1/4, or 0.25.
• ƒ(x) = 1/x, where x is [1,100], with 1,000 approximations. The exact result is the natural log of 100, or about 4.605170
• ƒ(x) = x, where x is [0,5000], with 5,000,000 approximations. The exact result is 12,500,000.
• ƒ(x) = x, where x is [0,6000], with 6,000,000 approximations. The exact result is 18,000,000.

源自於 https://rosettacode.org/wiki/Numerical_integration#C

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

double int_leftrect(double from, double to, double n, double (*func)())
{
   double h = (to-from)/n;
   double sum = 0.0, x;
   for(x=from; x <= (to-h); x += h)
      sum += func(x);
   return h*sum;
}

double int_rightrect(double from, double to, double n, double (*func)())
{
   double h = (to-from)/n;
   double sum = 0.0, x;
   for(x=from; x <= (to-h); x += h)
     sum += func(x+h);
   return h*sum;
}

double int_midrect(double from, double to, double n, double (*func)())
{
   double h = (to-from)/n;
   double sum = 0.0, x;
   for(x=from; x <= (to-h); x += h)
     sum += func(x+h/2.0);
   return h*sum;
}

double int_trapezium(double from, double to, double n, double (*func)())
{
   double h = (to - from) / n;
   double sum = func(from) + func(to);
   int i;
   for(i = 1;i < n;i++)
       sum += 2.0*func(from + i * h);
   return  h * sum / 2.0;
}

double int_simpson(double from, double to, double n, double (*func)())
{
   double h = (to - from) / n;
   double sum1 = 0.0;
   double sum2 = 0.0;
   int i;

   double x;

   for(i = 0;i < n;i++)
      sum1 += func(from + h * i + h / 2.0);

   for(i = 1;i < n;i++)
      sum2 += func(from + h * i);

   return h / 6.0 * (func(from) + func(to) + 4.0 * sum1 + 2.0 * sum2);
}


/* test */
double f3(double x)
{
  return x;
}

double f3a(double x)
{
  return x*x/2.0;
}

double f2(double x)
{
  return 1.0/x;
}

double f2a(double x)
{
  return log(x);
}

double f1(double x)
{
  return x*x*x;
}

double f1a(double x)
{
  return x*x*x*x/4.0;
}

typedef double (*pfunc)(double, double, double, double (*)());
typedef double (*rfunc)(double);

#define INTG(F,A,B) (F((B))-F((A)))

int main()
{
     int i, j;
     double ic;

     pfunc f[5] = {
       int_leftrect, int_rightrect,
       int_midrect,  int_trapezium,
       int_simpson
     };
     const char *names[5] = {
       "leftrect", "rightrect", "midrect",
       "trapezium", "simpson"
     };
     rfunc rf[] = { f1, f2, f3, f3 };
     rfunc If[] = { f1a, f2a, f3a, f3a };
     double ivals[] = {
       0.0, 1.0,
       1.0, 100.0,
       0.0, 5000.0,
       0.0, 6000.0
     };
     double approx[] = { 100.0, 1000.0, 5000000.0, 6000000.0 };

     for(j=0; j < (sizeof(rf) / sizeof(rfunc)); j++)
     {
       for(i=0; i < 5 ; i++)
       {
         ic = (*f[i])(ivals[2*j], ivals[2*j+1], approx[j], rf[j]);
         printf("%10s [ 0,1] num: %+lf, an: %lf\n",
                names[i], ic, INTG((*If[j]), ivals[2*j], ivals[2*j+1]));
       }
       printf("\n");
     }
}


輸出畫面
  leftrect [ 0,1] num: +0.235322, an: 0.250000
 rightrect [ 0,1] num: +0.245025, an: 0.250000
   midrect [ 0,1] num: +0.240137, an: 0.250000
 trapezium [ 0,1] num: +0.250025, an: 0.250000
   simpson [ 0,1] num: +0.250000, an: 0.250000

  leftrect [ 0,1] num: +4.654000, an: 4.605170
 rightrect [ 0,1] num: +4.555991, an: 4.605170
   midrect [ 0,1] num: +4.603772, an: 4.605170
 trapezium [ 0,1] num: +4.605986, an: 4.605170
   simpson [ 0,1] num: +4.605170, an: 4.605170

  leftrect [ 0,1] num: +12499992.500730, an: 12500000.000000
 rightrect [ 0,1] num: +12499997.500729, an: 12500000.000000
   midrect [ 0,1] num: +12499995.000729, an: 12500000.000000
 trapezium [ 0,1] num: +12500000.000000, an: 12500000.000000
   simpson [ 0,1] num: +12500000.000000, an: 12500000.000000

  leftrect [ 0,1] num: +17999991.001392, an: 18000000.000000
 rightrect [ 0,1] num: +17999997.001390, an: 18000000.000000
   midrect [ 0,1] num: +17999994.001391, an: 18000000.000000
 trapezium [ 0,1] num: +18000000.000000, an: 18000000.000000
   simpson [ 0,1] num: +18000000.000000, an: 18000000.000000