alex9ufo 聰明人求知心切: 12月 2018

2018年12月31日星期一

Simpson 辛普森積分 C語言與Python語言

Program for Simpson’s 1/3 Rule

源自於

https://www.geeksforgeeks.org/program-simpsons-13-rule/

In numerical analysis, Simpson’s 1/3 rule is a method for numerical approximation of definite integrals. Specifically, it is the following approximation:

$\int_{a}^{b} f(x) dx \approx \frac{(b-a)}{6} \bigg(f(a) + 4f \frac{(a+b)}{2} + f(b)\bigg)$

In Simpson’s 1/3 Rule, we use parabolas to approximate each part of the curve.We divide
the area into n equal segments of width Δx.
Simpson’s rule can be derived by approximating the integrand f (x) (in blue)
by the quadratic interpolant P(x) (in red).
graph

In order to integrate any function f(x) in the interval (a, b), follow the steps given below:

1.Select a value for n, which is the number of parts the interval is divided into.
2.Calculate the width, h = (b-a)/n
3.Calculate the values of x0 to xn as x0 = a, x1 = x0 + h, …..xn-1 = xn-2 + h, xn = b.
Consider y = f(x). Now find the values of y(y0 to yn) for the corresponding x(x0 to xn) values.
4.Substitute all the above found values in the Simpson’s Rule Formula to calculate the integral value.

Approximate value of the integral can be given by Simpson’s Rule:

$\int_{a}^{b} f(x) dx \approx \frac{h}{3} \bigg(f_0 + f_n + 4 * \sum_{i=1,3,5}^{n-1}f_i + 2* \sum_{i=2,4,6}^{n-2}f_i \bigg)$

Note : In this rule, n must be EVEN.

Application :
It is used when it is very difficult to solve the given integral mathematically.
This rule gives approximation easily without actually knowing the integration rules.

Example :

Evaluate logx dx within limit 4 to 5.2.

First we will divide interval into six equal 
parts as number of interval should be even.

x    :  4     4.2   4.4   4.6   4.8  5.0  5.2
logx :  1.38  1.43  1.48  1.52  1.56 1.60 1.64

Now we can calculate approximate value of integral
using above formula:
     = h/3[( 1.38 + 1.64) + 4 * (1.43 + 1.52 + 
                      1.60 ) +2 *(1.48 + 1.56)]
     = 1.84
Hence the approximation of above integral is 
1.827 using Simpson's 1/3 rule.

// CPP program for simpson's 1/3 rule 
#include <iostream> 
#include <math.h> 
using namespace std; 

// Function to calculate f(x) 
float func(float x) 
{ 
 return log(x); 
} 

// Function for approximate integral 
float simpsons_(float ll, float ul, int n) 
{ 
 // Calculating the value of h 
 float h = (ul - ll) / n; 

 // Array for storing value of x and f(x) 
 float x[10], fx[10]; 

 // Calculating values of x and f(x) 
 for (int i = 0; i <= n; i++) { 
  x[i] = ll + i * h; 
  fx[i] = func(x[i]); 
 } 

 // Calculating result 
 float res = 0; 
 for (int i = 0; i <= n; i++) { 
  if (i == 0 || i == n) 
   res += fx[i]; 
  else if (i % 2 != 0) 
   res += 4 * fx[i]; 
  else
   res += 2 * fx[i]; 
 } 
 res = res * (h / 3); 
 return res; 
} 

// Driver program 
int main() 
{ 
 float lower_limit = 4; // Lower limit 
 float upper_limit = 5.2; // Upper limit 
 int n = 6; // Number of interval 
 cout << simpsons_(lower_limit, upper_limit, n); 
 return 0; 
}

Output:

1.827847

Python 語言

# Python code for simpson's 1 / 3 rule 
import math 

# Function to calculate f(x) 
def func( x ): 
 return math.log(x) 

# Function for approximate integral 
def simpsons_( ll, ul, n ): 

 # Calculating the value of h 
 h = ( ul - ll )/n 

 # List for storing value of x and f(x) 
 x = list() 
 fx = list() 
 
 # Calcuting values of x and f(x) 
 i = 0
 while i<= n: 
  x.append(ll + i * h) 
  fx.append(func(x[i])) 
  i += 1

 # Calculating result 
 res = 0
 i = 0
 while i<= n: 
  if i == 0 or i == n: 
   res+= fx[i] 
  elif i % 2 != 0: 
   res+= 4 * fx[i] 
  else: 
   res+= 2 * fx[i] 
  i+= 1
 res = res * (h / 3) 
 return res 
 
# Driver code 
lower_limit = 4 # Lower limit 
upper_limit = 5.2 # Upper limit 
n = 6 # Number of interval 
print("%.6f"% simpsons_(lower_limit, upper_limit, n))

定積分之數值計算

源自於 http://www.stat.nuk.edu.tw/cbme/math/calculus/cal2/c6_4/bud.htm

定積分之數值計算

若不知 $f$ 之反導數, 則便無法利用微積分基本定理求定積分 $\int_a^b f(x)dx$ 之值。不過我們卻可求此定積分之近似值至其任意精確度。方法之一便是利用第二章定理 4.9, 以 Riemann 和來逼近。

在實際應用時, 可取 $P_1$ , $P_2,\cdots$ 為 $[a, b]$ 之一正規分割數列若 $P_n \linebreak =\{x_0,x_1,\cdots,x_n\}$ , 則可以

$\begin{displaymath} R(P_n)=\sum_{i=1}^n f(x_{i-1})\Delta x, \end{displaymath}$

或

$\begin{displaymath} R_1(P_n)=\sum_{i=1}^n f(x_i)\Delta x, \end{displaymath}$

當做第 $n$ 個 Riemann 和, 其中 $\Delta x=(b-a)/n$ 只要 $n$ 夠大, 可以 $R(P_n)$ 與 $R_1(P_n)$ 之算術平均來估計 $\int_a^b f(x)dx$ , 通常此為一更好的估計。而

$\begin{eqnarray*} &&\frac 1 2 (R(P_n)+R_1(P_n)) \\ &&=\frac {\Delta x} 2(\sum_{... ...\sum_{i=1}^{n-1}2f(x_i)+f(x_n))\mbox{\raisebox{-1.2mm}{\large }} \end{eqnarray*}$ 。

因此若 $P_n=\{x_0,x_1, \cdots,x_n\}$ 為 $[a, b]$ 之一正規分割, 則

$\begin{displaymath} \hspace*{0.8cm} \int_a^b f(x)dx \doteq \frac {b-a}{2n}(f(x_0... ...1)+\cdots+2f(x_{n-1})+f(x_n))\mbox{\raisebox{-1.2mm}{\large }} \end{displaymath}$ 。

(4.1)

上述積分公式便稱為梯形法 (trapezoidal rule)此名稱之由來如下

設 $f(x)\geq 0$ , $x\in [a, b]$ , 則如圖 4.1 可看出, (4.1) 式之右側為 $n$ 個高皆為 $\Delta x=(b-a)/n$ 之梯形的面積和

例 1.設 $f(x)=\sqrt {x^2+1}$ 利用梯形法求在 $f$ 之圖形下, 由 0 至 3 之面積。

其次來看梯形法之一推廣, 此法源自於 Simpson (1710-1761), 故稱之為 Simpson 法 (Simpson's rule), 為一較梯形法更精確的估計定積分的方法。

設 $f$ 為一在 $[a, b]$ 上連續之函數, 梯形法是以一線段來逼近 $f$ 之圖形, 而 Simpson 法是以一拋物線來逼近 $f$ 之圖形, 所以又稱拋物線法 (parabolic rule)。我們先看拋物線的一些性質給三個不共線之相異點 $(-\Delta x, y_0)$ , $(0, y_1)$ 及 $(\Delta x, y_2)$ , 恰有一拋物線

通過此三點, 其中係數 $a_0, a_1, a_2$ 滿足下述三方程式。

$\displaystyle y_0$	$\textstyle =$	$\displaystyle a_2(-\Delta x)^2+a_1(-\Delta x)+a_0,$	(4.2)
$\displaystyle y_1$	$\textstyle =$	$\displaystyle a_0,$	(4.3)
$\displaystyle y_2$	$\textstyle =$	$\displaystyle a_2 (\Delta x)^2+a_1\Delta x +a_0\mbox{\raisebox{-1.2mm}{\large }}$ 。	(4.4)

而在此拋物線下, 由 $-\Delta x$ 至之面積 (見圖 4.2) 為

$\begin{eqnarray*} \int_{-\Delta x}^{\Delta x}(a_2 x^2+a_1 x+a_0)dx &=& (\frac {a... ...lta x} 3(2a_2(\Delta x)^2+6a_0)\mbox{\raisebox{-1.2mm}{\large }} \end{eqnarray*}$ 。

由 (4.2)-(4.4), 可證明

$\begin{displaymath} y_0+4y_1+y_2=2a_2(\Delta x)^2+6a_0\mbox{\raisebox{-1.2mm}{\large }} \end{displaymath}$ 。

此即圖 4.2 中拋物線下的面積。

若將拋物線平移, 設有一拋物線 $y=a_0 x^2+a_1 x+a_0$ 通過 $(x_0, y_0)$ , $(x_1, y_1)$ 及 $(x_2, y_2)$ 等三點, 見圖 4.3, 其中 $\Delta x=x_2-x_1=x_1-x_0$ 則顯然拋物線下的面積仍相同, 即

。

設 $P_n=\{x_0,x_1, \cdots,x_n\}$ 為 $[a, b]$ 之一正規分割, 其中 $n$ 為一偶數, 且令 $\Delta x=(b-a)/n$ 。則

$\begin{eqnarray*} \int_a^b f(x)dx &=& \int_{x_0}^{x_2}f(x)dx+\int_{x_2}^{x_4}f(x... ... && +\int_{x_{n-2}}^{x_n}f(x)dx\mbox{\raisebox{-1.2mm}{\large }} \end{eqnarray*}$ 。

其中每一積分 $\int_{x_i}^{x_{i+2}}f(x)dx$ 皆可以一通過三點 $(x_i, f(x_i))$ , $(x_{i+1},\linebreak f(x_{i+1}))$ 及 $(x_{i+2}, f(x_{i+2}))$ 之拋物線下的面積來逼近又由前面已得的結果知:

$\begin{eqnarray*} \int_{x_0}^{x_2} f(x)dx &\doteq& \frac {\Delta x}{3}(f(x_0)+4f(x_1)+f(x_2)),\\ \end{eqnarray*}$

$\begin{eqnarray*} \int_{x_2}^{x_4} f(x)dx &\doteq& \frac {\Delta x}{3}(f(x_2)+4f... ...(f(x_{n-2})+4f(x_{n-1})+f(x_n))\mbox{\raisebox{-1.2mm}{\large }} \end{eqnarray*}$ 。

將這些式子左、右各分別相加, 即得

$\displaystyle \int_a^bf(x)dx$	$\textstyle \doteq$	$\displaystyle \frac {b-a}{3n}(f(x_0)+4f(x_1)+2f(x_2)+4f(x_3)+\cdots$
		$\displaystyle +2f(x_{n-2})+4f(x_{n-1})+f(x_n))\mbox{\raisebox{-1.2mm}{\large }}$ 。	(4.5)