alex9ufo 聰明人求知心切: Simpson’s 1/3rd Rule

2019年1月1日星期二

Simpson’s 1/3rd Rule

源自於 https://phoxis.org/2011/03/11/simpson-13-rd-rule/

A brief introduction to the Simpson’s 1/3rd rule and a uniform interval Composite Simpson’s 1/3rd Rule implementation.

Simpson’s 1/3rd Rule

The Simpson’s 1/3 $^{rd}$ rule is a numerical method to find the integral $\intop_{a}^{b}f(x)\, dx$ within some finite limits $a$ and $b$ . Simpson’s 1/3rd rule approximates $f(x)$ with a polynomial of degree two $P(x)$ , ie. a parabola between the two limits $a$ and $b$ , and then finds the integral of that bounded parabola, and is used to represent the approximate integral $\intop_{a}^{b}f(x)\, dx$ . The integral of the approximated function is the area under the parabola bounded by the points $a$ and $b$ and by the positive side of the x axis. The quadratic function has three points common to the function $f(x)$ , which are as follows: The end points of the approximate quadratic function $P(x)$ is the same as the function $f(x)$ at $a$ , $b$ . And $P(x)$ takes the same value of the function $f(x)$ at point $m = (a+b)/2$ . Thus three points are fixed each in equal interval $a>m>b$ and a parabola is drawn through these three points $f(a),f(m),f(b)$ and the area under the parabola through these points bounded by $a$ and $b$ and the positive side of the X axis is found, which is used as the approximated integral value. The below iterative formula an be used to find the integral of a function $f(x)$ using Simpson’s 1/3rd rule.

$I_{S1}=\frac{h}{3}(y_{0}+4y_{1}+y_{2})\mbox{ where }y_{i}=f(x_{i})$

The above formula approximates $f(x)$ with only one parabola. To represent the function more accurately with with method, the function is divided in many intervals say $n$ having same width say $h$ , such that $h=\frac{(b-a)}{n}$ and then a parabola is fit within each of the $n$ intervals to approximate the integrals of the function within those intervals. The evaluation of the integral within each such interval using the above single interval iterative formula is used to find the overall evaluation of the integral within the limits $a$ and $b$ . This is known as the Simpson’s 1/3rd Composite Rule. The iterative method for the Simpson’s 1/3rdComposite Rule can be found by repetitively applying the single interval formula in each intervals, and is found to be:

$I_{S}=\frac{h}{3}[(y_{0}+y_{n})+2\sum_{i=even}y_{i}+4\sum_{j=odd}y_{j}]$
$\mbox { where }y_{i}=f(x_{i}),y_{j}=f(x_{j})\mbox{ and }(0 < i < n) , (0 < j < n)$

The Process

The Simpson’s $1/3{rd}$ Composite Rule function will calculate the first and the last ordinate separately and also calculate and sum up the odd and the even ordinates separately for better clarity, and at last the Simpson’s $1/3{rd}$ Composite Rule iterative formula would be used directly.
For demonstration purpose some sample functions are used:

$f(x)=\frac{1}{x}+5+10x^{2}$
$f(x)=\frac{sin(2x)}{(1+x)^{5}}$

FUNCTION simpson_13(f, a, b, n)

  Find width of each interval h = (b-a)/n 

  Find y0 = f(a), yn = f(b)

  Initialize x = a+h, i = 1, sum = 0

  Initialize even = 0, odd = 0

  /* Calculate the sum of the even and odd ordinates from (x + h) to (x + (n-1)h) */

  WHILE (i < n)

   DO

    IF i is even

     THEN

      even = even + f (x)

    ELSE IF i is odd

     THEN

      odd = odd + f(x)

    ENDIF

    x = x + h

  ENDWHILE

  /* Apply Simpson's 1/3rd Composite Formula to find the integral */

  Is = (h/3* {y0 + yn + 2 * even + 4 * odd})

ENDFUNCTION

Sourcecode

#include <stdio.h>

#include <math.h>

double f1 (double x);

double f2 (double x);

double simpson_13 (double (*f) (double x), double a, double b, int n);

int

main (void)

{

  double a, b, n;

  double Is;

  printf ("\nEnter a,b,n: ");

  scanf ("%lf %lf %lf", &a, &b, &n);

  printf ("\nf(x) = sin (2x) / (1+x)^5");

  /* Show integral computed with simpson's 1/3 rd rule */

  Is = simpson_13 (f1, a, b, n);

  printf ("\n\tI_Simpson_13rd (f(x), %g, %g, %g) = %g", a, b, n, Is);

  printf ("\n");

  printf ("\nf(x) = (1/x) + 5 + 10x^2;");

  /* Show integral computed with simpson's 1/3 rd rule */

  Is = simpson_13 (f2, a, b, n);

  printf ("\n\tI_Simpson_13rd (f(x), %g, %g, %g) = %g", a, b, n, Is);

  printf ("\n");

  return 0;

}

/* Takes a function pointer which should represent the single variable

 * function of which the integral is to be found. The function should

 * receive a double and return a double. The limits a and b and the

 * number of intervals n is provided. The value returned is the integral

 * of the passed function is the argument

 */

double

simpson_13 (double (*f) (double x), double a, double b, int n)

{

  double h;

  double y = 0, x, even = 0, odd = 0, y0, yn;

  int i;

  h = (b - a) / n;

  y0 = f (a);

  yn = f (b);

  for (i = 1, x = a + h; i < n; x = x + h, i++)

    {

      if (i % 2 == 0)

        even = even + f (x);

      else

        odd = odd + f (x);

    }

  y = (h / 3) * ((y0 + yn) + 2 * even + 4 * odd);

  return y;

}

double

f1 (double x)

{

  return sin (2 * x) / pow ((1 + x), 5);

}

double

f2 (double x)

{

  return (1 / x) + 5 + 10 * x * x;

}

Description

The two given function defined in the previous section is implemented with the functions f1 () and f2 (). These two functions accept a double type data, and returns double.

double f1 (double x) : This function implements the equation sin (2 * x) / pow ((1 + x), 5) and returns the evaluated value.
double f2 (double x) : This function implements the equation (1 / x) + 5 + 10 * x * x and returns the evaluated value.
double simpson_13 (double (*f) (double x), double a, double b, int n) : This takes a pointer to a function pointer which accepts one double type data and returns a double type data. The next argument is double a, which is interpreted as the lower limit of the integral, the next argument double b is considered as the upper limit of the integral. int n is the number of intervals to be taken. After finding the width of each interval, this first calculates the first and the last ordinate y0 = f(a) and y1 = f(b). Then it enters the for loop and computes the $(n-2)$ values of the passed function starting from $(x+h)$ up to $(x+(n-1)h$ ). The values of the function values for the odd interval number are summed in a separate variable odd, and the values for the even interval numbers are summed in another variable even. At the end of the interval the composite iterative Simpson’s formula is directly used to find the integral, and the value is returned:

(h / 3) * ((y0 + yn) + 2 * even + 4 * odd);
int main (void) : The main function prompts the user to enter the upper and the lower limits and the number of intervals to be taken and calls the functions with the proper parameters. The simpson_13() is called as Is = simpson_13 (f1, a, b, n);which computes the integral of the function f1 within upper and lower limit a and b with the Simpson’s 1/3 $^{rd}$ rule.

Error

When $f^{(4)}$ is continuous and $M$ is the upper bound for the values
of $\left|f^{(4)}\right|$ on $[a,b]$ , then Simpson’s 1/3rd is found to be: $\left|E_{S}\right|\leq\frac{(b-a)}{180}h^{4}M$

The magnitude of the Simpson’s rule error decreases as the the step size is decreased as the function is better approximated with the second degree equation.

Comments

The Simpson’s 1/3rd evaluation function is implemented to be independent of the function of which the integral is to be found, by allowing passing of the function to be integrated as a parameter. This implementation could be incorporated into an algebraic equation parser, which can be used to evaluate some user entered functions to be evaluated like for example:
ans = simpson_13 (parse_expression (user_expression, x), a, b, n) ; the ‘user_expression’ is being considered a single variable function with the variable replaced with the value of ‘x’ .

Output

Some sample outputs of the above code is shown below

Run 1: input = lower limit a = 1, upper limit b = 2, intervals n = 10

Enter a,b,n: 1 2 10

f(x) = sin (2x) / (1+x)^5

    I_Simpson_13rd (f(x), 1, 2, 10) = 0.00505831

f(x) = (1/x) + 5 + 10x^2;

    I_Simpson_13rd (f(x), 1, 2, 10) = 29.0265

Run 2: input = lower limit a = 5.3, upper limit b = 10.23, intervals n = 100

Enter a,b,n: 5.3 10.23 100

f(x) = sin (2x) / (1+x)^5

    I_Simpson_13rd (f(x), 5.3, 10.23, 100) = -3.20359e-05

f(x) = (1/x) + 5 + 10x^2;

    I_Simpson_13rd (f(x), 5.3, 10.23, 100) = 3097.71

alex9ufo 聰明人求知心切

2019年1月1日星期二

Simpson’s 1/3rd Rule

Simpson’s 1/3rd Rule

Simpson’s 1/3rd Rule

The Process

Sourcecode

Description

Error

Comments

Output

Links

References

沒有留言:

張貼留言

WOKWI LED + MQTT Node-Red SQLite

2019年1月1日 星期二

Simpson’s 1/3rd Rule

Simpson’s 1/3rd Rule

Simpson’s 1/3rd Rule

The Process

Sourcecode

Description

Error

Comments

Output

Links

References

沒有留言:

張貼留言

WOKWI LED + MQTT Node-Red SQLite

2019年1月1日星期二