6. 積分

課程筆記

積分 (Integration)

Newton-Cotes 積分法

• 最常用的數值積分法。
• 原理：
  2. 其中， 為 次多項式：
  3. 可改寫為
     = + + + + +
     其中， = ， 為數據點間距。
• 誤差：
• 多重梯形法：
  • 誤差：，其中 為 的平均值。
• 討論：
  • 多區段梯形法能提供足夠的準確度。
  • 高準確度要求下，多區段梯形法需更多計算。
  • 捨入誤差影響也限制計算積分值得能力。

辛普森法則 (Simpson’s Rule)

• 辛普森 1/3 法則：積分值為位於連接三個點的拋物線下方的區域面積。
  • 積分：
  • 誤差：
    • ，其中
  • 原理：
    2. 以 、 代替 、
    3. ，
    4. 得：
  • 多重積分（區段數必須為偶數）：
    • 誤差：
• 辛普森 3/8 法則：積分值為位於連接四個點的三階多項式下方的區域面積。
  • 積分：
  • 誤差：
    • ，其中
• 辛普森 1/3 法則用三點可達三階準確度，而 3/8 法則用四點才能達到三階準確度。因此，1/3 法則較為常用。但是，3/8 法則可以用用在積分區間為奇數時。因此，一般當積分區間為奇數時，會先使用 1/3 法則，再使用 3/8 法則處理最後三個區間。
  • 若區段數為偶數，每兩區段套用 1/3 法則。
  • 若區段數為奇數，每兩區段套用 1/3 法則，最後三區段套用 3/8 法則。

不等距離區間上的積分 (Integration With Unequal Intervals)

策略：

• 使用梯形法計算不等距離區間的積分值。
• 等間距時使用辛普森 1/3、3/8 法則。

Example

在 與 之間的積分值。

Solution-Simpson’s 1/3 Rule

1. ；；
3. 誤差：；
4. 誤差估計值：

Solution-Simpson’s 1/3 多重積分

1. ；；；；
3. 誤差：；
4. 誤差估計值：

Solution-Simpson’s 3/8 多重積分

1. ；；；；；
2. 前二個區段使用 1/3 法則：
3. 最後三個區段使用 3/8 法則：
5. 誤差：；

重積分 (Multiple Integration)

• 二維函數平均值：

Example

，在 與 之間的積分值。
其中，已知九點結果如下：

Solution-梯形法

1. 梯形法公式：
2. 列處理：
   1. ；；，因此
   2. 第一列：
   3. 第二列：
   4. 第三列：
3. 行處理 (整體)：
   1. ；；，因此
   2. 積分值：
4. 平均值：

Solution-Simpson’s 1/3 Rule

1. 由於 ，因此使用 Simpson’s 1/3 Rule。
2. 辛普森 1/3 法則：
3. 列處理：
   1. ；；，因此
   2. 第一列：
   3. 第二列：
   4. 第三列：
4. 行處理 (整體)：
   1. ；；，因此
   2. 積分值：
5. 平均值：

Homework

題目

試使用梯形法、Simposon’s 1/3 法則計算：

1. 當 時，計算積分值。
2. 當 時，計算積分值。
3. 當 時，計算積分值。

解題流程

1. 首先，求得函式的真值為 8/3，詳如式 (1)：

2. 於本次作業中，由於 n 皆為偶數倍，因此可全部使用 Simposon’s 1/3 法則評估，不須對奇數區間部分使用 Simposon’s 3/8 方法。

Code

針對本作業撰寫 Golang 代碼如下所示。代碼分為四個核心部分，分別為 defPoints()、calcTrapezoidal()、calcSimpson()、Homework()。為利於說明，從 Homework() 方法開始，該方法定義了函式 (1) 為 f，以及其積分區間（x0、x1、y0、y1）與真值（trueValue）。並根據題目要求，分別使用 n=10、20、50 進行計算。
在計算前，透過 defPoints() 方法根據積分區間、函式以及區間數量建立 (n+1) * (n+1) * 3 形式的三維矩陣，分別存放 x, y, f(x, y)。以作為梯形法與辛普森方法的數據輸入（points）。
而 calcTrapezoidal()、calcSimpson() 分別為梯形法與辛普森法的實現。兩者開頭皆相同，皆為計算 h (Δx)、k (Δy)，從輸入數據讓方法了解輸入數據的間隔。接著進行計算。梯形法部分計算 sum = f0 + 2Σfi + fn 再透過 I = h * k * sum / 4 求得積分值。而辛普森法則是計算 sum = f0 + 4Σfi奇 + 2Σfi偶 + fn （先建立權重清單 simpsonWeight，再計算加權）再透過 I = h * k * sum / 9 求得積分值。最後，計算真值與積分值之差值與相對誤差。為了觀察到誤差分布，本作業將誤差顯示至小數點下 16 位。

package main
 
import (
    "fmt"
)
 
func defPoints(x0 float64, x1 float64, y0 float64, y1 float64, f func(float64, float64) float64, n int) [][][]float64 {
    h := (x1 - x0) / float64(n)
    k := (y1 - y0) / float64(n)
 
    points := make([][][]float64, n+1)
    for i := 0; i <= n; i++ {
        points[i] = make([][]float64, n+1)
        for j := 0; j <= n; j++ {
            xVal := x0 + float64(i)*h
            yVal := y0 + float64(j)*k
            points[i][j] = []float64{xVal, yVal, f(xVal, yVal)}
        }
    }
 
    return points
}
 
func calcTrapezoidal(points [][][]float64, trueValue float64) {
    n := len(points) - 1
 
    h := points[1][0][0] - points[0][0][0]
    k := points[0][1][1] - points[0][0][1]
 
    var sum float64
    for i := 0; i <= n; i++ {
        for j := 0; j <= n; j++ {
            weight := 1.0
            if i > 0 && i < n {weight *= 2}
            if j > 0 && j < n {weight *= 2}
            sum += weight * points[i][j][2]
        }
    }
    integral := (h * k / 4.0) * sum
    err := trueValue - integral
    errt := (trueValue - integral) / trueValue * 100
    fmt.Printf("Trapezoidal: %.16f, Error: % .16f, Error rate: % .16f\n", integral, err, errt)
}
 
func calcSimpson(points [][][]float64, trueValue float64) {
    n := len(points) - 1
 
    h := points[1][0][0] - points[0][0][0]
    k := points[0][1][1] - points[0][0][1]
 
    simpsonWeight := func(idx, n int) float64 {
        if idx == 0 || idx == n {
            return 1.0
        } else if idx%2 == 1 {
            return 4.0
        } else {
            return 2.0
        }
    }
 
    var sum float64
    for i := 0; i <= n; i++ {
        for j := 0; j <= n; j++ {
            wx := simpsonWeight(i, n)
            wy := simpsonWeight(j, n)
            sum += wx * wy * points[i][j][2]
        }
    }
 
    integral := (h * k / 9.0) * sum
    err := trueValue - integral
    errt := (trueValue - integral) / trueValue * 100
    fmt.Printf("Trapezoidal: %.16f, Error: % .16f, Error rate: % .16f\n", integral, err, errt)
}
 
func Homework() {
    f := func(x float64, y float64) float64 {
        return x * x + y * y
    }
    x0 := -1.0
    x1 :=  1.0
    y0 := -1.0
    y1 :=  1.0
 
    trueValue := 8.0 / 3.0
 
    // n=10
    fmt.Printf(  "=== n = 10 ===\n")
    points := defPoints(x0, x1, y0, y1, f, 10)
    calcTrapezoidal(points, trueValue)
    calcSimpson(points, trueValue)
 
    // n=20
    fmt.Printf("\n=== n = 20 ===\n")
    points = defPoints(x0, x1, y0, y1, f, 20)
    calcTrapezoidal(points, trueValue)
    calcSimpson(points, trueValue)
    
    // n=50
    fmt.Printf("\n=== n = 50 ===\n")
    points = defPoints(x0, x1, y0, y1, f, 50)
    calcTrapezoidal(points, trueValue)
    calcSimpson(points, trueValue)
}
 
func main() {
    Homework()
}

---

< Last

Next >