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

References

Calculus and Analytic Geometry: Thomas Finney
Images from : Wikipedia

辛普森積分法則

Simpson's rule 辛普森積分法則

From Wikipedia, the free encyclopedia

Jump to navigation Jump to search

\int _{a}^{b}f(x)\,dx\approx {\tfrac {\Delta x}{3}}\left(f(x_{0})+4f(x_{1})+2f(x_{2})+4f(x_{3})+2f(x_{4})+\cdots +4f(x_{n-1})+f(x_{n})\right)

,In numerical analysis, Simpson's rule is a method for numerical integration, the numerical approximation of definite integrals. Specifically, it is the following approximation for

n

equally spaced subdivisions (where

n

is even): (General Form)

where

\Delta x={\frac {b-a}{n}}

and

x_{i}=a+i\Delta x

//Simpson's 1/3rd Rule for Evaluation of Definite Integrals
#include<iostream>
#include<cmath>
using namespace std;
double f(double x)
{
double a=1/(1+x*x); //write the function whose definite integral is to be calcuated here
return a;
}
int main()
{ cout.precision(4); //set the precision
cout.setf(ios::fixed);
int n,i; //n is for subintervals and i is for loop
double a,b,c,h,sum=0,integral;
cout<<"\nEnter the limits of integration,\n\nInitial limit,a= ";
cin>>a;
cout<<"\nFinal limit, b="; //get the limits of integration
cin>>b;
cout<<"\nEnter the no. of subintervals(IT SHOULD BE EVEN), \nn="; //get the no. of subintervals
cin>>n;
double x[n+1],y[n+1];
h=(b-a)/n; //get the width of the subintervals
for (i=0;i<n+1;i++)
{ //loop to evaluate x0,...xn and y0,...yn
x[i]=a+i*h; //and store them in arrays
y[i]=f(x[i]);
}
for (i=1;i<n;i+=2)
{
sum=sum+4.0*y[i]; //loop to evaluate 4*(y1+y3+y5+...+yn-1)
}
for (i=2;i<n-1;i+=2)
{
sum=sum+2.0*y[i]; /*loop to evaluate 4*(y1+y3+y5+...+yn-1)+
2*(y2+y4+y6+...+yn-2)*/
}
integral=h/3.0*(y[0]+y[n]+sum); //h/3*[y0+yn+4*(y1+y3+y5+...+yn-1)+2*(y2+y4+y6+...+yn-2)]
cout<<"\nThe definite integral is "<<integral<<"\n"<<endl;
return 0;
}

======================

輸出畫面

======================

Enter the limits of integration,

Initial limit,a= 0

Final limit, b=6

Enter the no. of subintervals(IT SHOULD BE EVEN),

n=10

The definite integral is 1.4002

數值積分梯形法辛普森1/3法則及3/8法則

http://keejko.blogspot.com/2014/08/blog-post.html

數值積分

許多函數的積分結果並沒有解析解，此時就要利用數值積分算出近似值，常用的方法有

矩形法（維基百科　中文　en）
梯形法（維基百科　中文　en）
辛普森1/3法則及3/8法則（維基百科　中文　en）

為了完成這些計算，我們可以自己寫程式，也可以利用LibreOffice Calc、Microsoft Excel這類的工作表軟體，甚至還可以利用Wolfram Alpha網站(http://www.wolframalpha.com/)查詢答案。以運算速度以及可以分割的區間數量而言，自己寫程式一定是最好的。

以下是我用LibreOffice Calc製作的檔案：數值積分.ods

另外我也用c語言寫了幾個小程式，主要是從維基百科上找到的程式碼改寫的，測試的環境為Windows 8.1 64位元版本、Dev-C++ 5.7.1版，以下是原始碼：

/* 數值積分：梯形法

由維基百科上的程式碼改寫而成

來源： http://zh.wikipedia.org/wiki/矩形法

改作者：王一哲 Yi-Zhe Wang

日期：Aug. 21, 2014 */

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

// 定義函數f(x)，也可以回傳其它數學函數

double f(double x){

return sin(x);

}

// 定義運算式

double trapezoidal_integrate(double a, double b, int subintervals){

double result;

double interval;

int i;

interval=(b-a)/subintervals;

result=0;

//首項、末項係數為1，其餘係數為2

for(i=1;i<=subintervals;i++){

if((i==1)||(i==subintervals)){

result+=f(a+interval*(i-1));

} else {

result+=2*f(a+interval*(i-1));

}

result*=interval/2;

return result;

}

//主程式

int main(void){

double integral;

integral=trapezoidal_integrate(0,2,500000);

printf("Integral: %.15f \n",integral);

system("pause");

return 0;

}

/* 數值積分：辛普森1/3法則

由維基百科上的程式碼改寫而成

來源： http://zh.wikipedia.org/wiki/矩形法

改作者：王一哲 Yi-Zhe Wang

日期：Aug. 21, 2014 */

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

// 定義函數f(x)，也可以回傳其它數學函數

double f(double x){

return sin(x);

}

// 定義運算式

double simpson1_3_integrate(double a, double b, int subintervals){

double result;

double interval;

int i;

interval=(b-a)/subintervals;

result=0;

//首項、末項係數為1，其餘每2項1組，第1項係數為4，第2項係數為2

for(i=0;i<=subintervals;i++){

if((i==0)||(i==subintervals)){

result+=f(a+interval*i);

} else if((i%2)==1){

result+=4*f(a+interval*i);

} else if((i%2)==0){

result+=2*f(a+interval*i);

}

result*=interval/3;

return result;

}

//主程式，分割的區間數只能是偶數

int main(void){

double integral;

integral=simpson1_3_integrate(0,2,500000);

printf("Integral: %.15f \n",integral);

system("pause");

return 0;

}

/* 數值積分：辛普森3/8法則

由維基百科上的程式碼改寫而成

來源： http://zh.wikipedia.org/wiki/矩形法

改作者：王一哲 Yi-Zhe Wang

日期：Aug. 21, 2014 */

#include <stdio.h>

#include <math.h>

#include <stdlib.h>

// 定義函數f(x)，也可以回傳其它數學函數

double f(double x){

return sin(x);

}

// 定義運算式

double simpson3_8_integrate(double a, double b, int subintervals){

double result;

double interval;

int i;

interval=(b-a)/subintervals;

result=0;

//首項、末項係數為1，其餘每3項1組，1、2項係數為3，第3項係數為2

for(i=0;i<=subintervals;i++){

if((i==0)||(i==subintervals)){

result+=f(a+interval*i);

} else if(((i%3)==1)||((i%3)==2)){

result+=3*f(a+interval*i);

} else if((i%3)==0){

result+=2*f(a+interval*i);

}

result*=interval*3/8;

return result;

}

//主程式，分割的區間數只能是3的倍數
int main(void){
double integral;
integral=simpson3_8_integrate(0,2,500001);

printf("Integral: %.15f \n",integral);

system("pause");

return 0;

}

如果要利用Wolfram Alpha網站查詢答案，就在方框內輸入

Integrate[f(x),{x,a,b}]

上式中的f(x)請換成要積分的函數，a換成積分下限，b換成積分下限，最後按下Enter即可。但如果想要看到進一步的資料，就必須付費註冊才行。

習題 EX4-2 梯行法雙重積分---C語言

輸入資料
n=10 a=0 b=5.0
0 0
0.5 4.67
1.0 7.34
1.5 8.86
2.0 9.73
2.5 10.22
3.0 10.51
3.5 10.67
4.0 10.76
4.5 10.81
5.0 10.81

程式
/* p4-2.c based on Simpson's Rule to compute
* definite integral with domain [a,b] and
* n even-grid. n must be even.
*/
#include <stdio.h>
#include <math.h>
void main()
{
int i,m,n;
double a,b,h,x[50],y[50],sum1=0.0,sum2=0.0,sn;
scanf("n=%d a=%lf b=%lf",&n,&a,&b);
for(i=0;i<=n;i++)
scanf("%lf %lf",&x[i],&y[i]);
m=n/2;
h=(b-a)/n;
for(i=1;i<=2*m-1;i++)
{
if(i%2==0)
sum2=sum2+y[i];
else
sum1=sum1+y[i];

sn=(h/3.0)*(y[0]+y[n]+2.0*sum2+4.0*sum1);
printf("S%d=%lf\n",n,sn);

}
//sn=(h/3.0)*(y[0]+y[n]+2.0*sum2+4.0*sum1);
//printf("S%d=%lf\n",n,sn);
return;
}
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

輸出結果

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
S10=4.915000
S10=7.361667
S10=13.268333
S10=16.511667
S10=23.325000
S10=26.828333
S10=33.941667
S10=37.528333
S10=44.735000

習題EX4-1 利用梯行法解

習題EX4-1 利用梯行法解

n=11 a=1 b=3.2
x y
======
1.0 4.9
1.2 5.4
1.4 5.8
1.6 6.2
1.8 6.7
2.0 7.0
2.2 7.3
2.4 7.5
2.6 8.0
2.8 8.2
3.0 8.3
3.2 8.3

/* p4-1.c based on Trapezoidol Rule is
* used for computing definite integral with
* domain [a,b] with n even-grids.
*/
#include <stdio.h>
#include <math.h>
#define PI 3.141593
void main()
{
int n,i;
double a,b,x[50],y[50],tn,h ,sum=0.0;
scanf("n=%d a=%lf b=%lf",&n,&a,&b);
for(i=1;i<=n+1;i++)
scanf("%lf %lf",&x[i],&y[i]);

h=(b-a)/n;
for(i=2;i<=n;i++)
{
sum=sum+y[i];
tn=(h/2.0)*(y[1]+y[n+1]+2.0*sum);
printf("T%d=%10.6lf\n",n,tn);
}

//tn=(h/2.0)*(y[1]+y[n+1]+2.0*sum);
//printf("T%d=%10.6lf\n",n,tn);
return;
}

T11= 2.400000
T11= 3.560000
T11= 4.800000
T11= 6.140000
T11= 7.540000
T11= 9.000000
T11= 10.500000
T11= 12.100000
T11= 13.740000
T11= 15.400000

訂閱：文章 (Atom)

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

辛普森積分法則

Simpson's rule 辛普森積分法則

數值積分梯形法辛普森1/3法則及3/8法則

數值積分

習題 EX4-2 梯行法雙重積分---C語言

習題EX4-1 利用梯行法解

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

辛普森積分法則

Simpson's rule 辛普森積分法則

數值積分 梯形法 辛普森1/3法則及3/8法則

數值積分

習題 EX4-2 梯行法雙重積分---C語言

習題EX4-1 利用梯行法解

2019年1月1日星期二

數值積分梯形法辛普森1/3法則及3/8法則