LoginSignup
2
1

More than 5 years have passed since last update.

@sort

数値積分を意外と早くないGo言語で書いてみた

Last updated at Posted at 2018-11-22

数値積分とは

積分区間[a,b]を等間隔に分割し、各区間での関数を多項式で近似して定積分を求める数値計算法の一つである。
以下に、定積分の定義を示す。
$$ \int_a^bf(x)dx \approx \lim_{\Delta x \to 0}\sum_{k=a}^bf(t_k)\Delta x $$
この近似式を使うことで近似的に積分を表現する。

数値積分の方法とソースコード

矩形法

矩形法とは

矩形法は領域を細長い長方形の短冊で刻んで足し、面積を求める。今回の場合は、短冊の左隅の値を用いるLeft-Ruleで実装する。

ソースコード

LeftRule.go
package main

import (
    "fmt"
    "math"
)

const (
    A = 0.0
    B = 1.0
    N = 262144 //max partition
    //TRUE = 3.0 /* 3*x^2+2*x+1 */
    //TRUE = 0.63661977236758134307663505349006 /* sin(x*pi/2) */
    //TRUE = 2.0 /* 2*x+1 */
    TRUE = 0.63661977236758134307553505349006 /* sin(x*pi) */
)

func integration(i float64) {
    dx :=0.0 
    s := 0.0
    err := 0.0
    dx = ( B - A )/i
    for x := A; x < B; x+=dx {
        //s += math.Sin( x*math.Pi/2 )*dx
        //s += ( 3*math.Pow( x, 2.0 ) + 2*x + 1 )*dx
        //s += ( 2*x + 1 )*dx
        s += math.Sin( x*math.Pi )*dx
    }
    err = math.Abs(( s - TRUE )/TRUE)
    fmt.Printf("division number = %.f\n",i)
    fmt.Printf("integral value = %.22f\n",s)
    fmt.Printf("relative error = %e\n\n",err)
}

func main() {
    for i := 2.0; i <= N; i *= 2 {
        integration(i)
    }
}

台形則

台形則とは

台形則とは定積分の値を細かく区切った台形の面積の和で近似する方法である。

ソースコード

TrapezoidMethod.go
package main

import (
    "fmt"
    "math"
)

const (
    A = 0.0
    B = 1.0
    N = 262144 //max partition
    //TRUE = 3.0 /* 3*x^2+2*x+1 */
    //TRUE = 0.63661977236758134307663505349006 /* sin(x*pi/2) */
    //TRUE = 2.0 /* 2*x+1 */
    TRUE = 0.63661977236758134307553505349006 /* sin(x*pi) */
)

func f(x float64) float64 {
    //return 3*math.Pow( x, 2.0 ) + 2*x + 1
    //return math.Sin(x*math.Pi/2)
    //return 2*x + 1
    return math.Sin(x*math.Pi)
}

func integration(i float64) {
    x := 0.0
    x1 := 0.0
    dx :=0.0 
    s := 0.0 
    s1 := 0.0
    err := 0.0
    dx = (B - A)/i
    for j := 0.0; j < i; j++ {
        x = A + j*dx
        x1 = x + dx
        s1 = ( ( ( f( x ) ) + ( f( x1 ) ) )/2 )*dx
        s += s1
    }
    err = math.Abs(( s - TRUE )/TRUE)
    fmt.Printf("division number = %.f\n",i)
    fmt.Printf("integral value = %.22f\n",s)
    fmt.Printf("relative error = %e\n\n",err)
}

func main() {
    for i := 2.0; i <= N; i *= 2 {
        integration(i)
    }
}

シンプソン則

シンプソン則とは

シンプソン則は関数を2次のラグランジュ補間によって定積分の値を求める方法である。

ソースコード

Simpson.go
package main

import(
    "fmt"
    "math"
)

const (
    A = 0.0
    B = 1.0
    N = 262144 //max partition
    TRUE = 0.63661977236758134307663505349006
)

func f(x float64) float64 {
    return math.Sin((x*math.Pi)/2)
}

func integration(i float64) {
    s := 0.0
    dx := 0.0
    err := 0.0
    dx = ( B - A )/i
    for x := A ; x < B ; x += dx*2 {
        s += ( f(x) + 4 * f(x + dx) + f(x + dx*2) )
    }
    s *= (B - A) * (dx / 3)
    err = math.Abs(( s - TRUE )/TRUE)
    fmt.Printf("division number = %.f\n",i)
    fmt.Printf("integral value = %.22f\n",s)
    fmt.Printf("relative error = %e\n\n",err)
}

func main() {
    for i := 2.0; i <= N; i *= 2 {
        integration(i)
    }
}

ロンバーグ則

ロンバーグ則とは

ロンバーグ則とは、定積分を合成台形公式で近似し、さらにRichardsonの補外法で近似をより正確にする方法である。

ソースコード

Romberg.go
package main

import (
    "fmt"
    "math"
    "./cast"
)

const (
    A = 0.0
    B = 1.0
    N = 18 //max partition 2**N
    TRUE = 0.63661977236758134307663505349006
)

var (
    s [N+1][N+1] float64
)

func f(x float64) float64 {
    return math.Sin( x*math.Pi/2 )
}

func p(x float64) float64 {
    return math.Pow( 4.0, x )
}

func trapezoidintegration(i int) {
    x := 0.0
    x1 := 0.0
    dx :=0.0 
    s1 := 0.0
    n := math.Pow(2.0, cast.FLOAT(i))
    dx = (B - A)/n
    for j := 0.0; j < n; j++ {
        x = A + j*dx
        x1 = x + dx
        s1 = ( ( f(x) + f(x1) )/2 )*dx
        s[i][1] += s1
    }
}

func integration() {
    for i := 2; i <= N; i++ {
        s[i][2] = ( 4*s[i][1] - s[i-1][1] )/3
        for j := 3; j <= i; j++ {
            s[i][j] = ( p( cast.FLOAT(j) )*s[i][j-1] - s[i-1][j-1] )/ ( p( cast.FLOAT(j) ) - 1 )
        }
    }
} 

func print() {
    var err float64
    for i := 1; i <= N; i++ {
        for j := 1; j <= i; j++ {
            fmt.Printf("%.22f ", s[i][j])
        }
        fmt.Printf("\n")
    }
    fmt.Printf("\n\n")
    for i := 1; i <= N; i++ {
        err = math.Abs(( s[i][i] - TRUE )/TRUE)
        fmt.Printf("division number = %d ",i)
        fmt.Printf("integral value = %.22f\n",s[i][i])
        fmt.Printf("relative error = %e\n",err)
    }
}

func main() {
    for i := 1; i <= N; i++ {
        trapezoidintegration(i)
    }
    integration()
    print()
}

実行結果

以下に書く方法で$\sin(x\pi/2)$の範囲[0,1]での定積分を求めた結果と真値との相対誤差を示す。

矩形法

LeftRule.exe
division number = 2
integral value = 0.3535533905932737308575
relative error = 4.446396e-01

division number = 4
integral value = 0.5034174365157310093721
relative error = 2.092337e-01

division number = 8
integral value = 0.5720731492255537453673
relative error = 1.013896e-01

division number = 16
integral value = 0.6048583632808495202937
relative error = 4.989070e-02

division number = 32
integral value = 0.6208669355013015644928
relative error = 2.474450e-02

division number = 64
integral value = 0.6287753141136420698132
relative error = 1.232205e-02

division number = 128
integral value = 0.6327055328642574583142
relative error = 6.148473e-03

division number = 256
integral value = 0.6346646499955104214763
relative error = 3.071099e-03

division number = 512
integral value = 0.6356427105247994813197
relative error = 1.534765e-03

division number = 1024
integral value = 0.6361313662819000347426
relative error = 7.671865e-04

division number = 2048
integral value = 0.6363756005336611565326
relative error = 3.835442e-04

division number = 4096
integral value = 0.6364976942528528525145
relative error = 1.917599e-04

division number = 8192
integral value = 0.6365587352607717797071
relative error = 9.587686e-05

division number = 16384
integral value = 0.6365892543018182303527
relative error = 4.793767e-05

division number = 32768
integral value = 0.6366045134566075125448
relative error = 2.396864e-05

division number = 65536
integral value = 0.6366121429425779520983
relative error = 1.198427e-05

division number = 131072
integral value = 0.6366159576626967409041
relative error = 5.992124e-06

division number = 262144
integral value = 0.6366178650170443154011
relative error = 2.996059e-06

台形則

TrapezoidMethod.exe
division number = 2
integral value = 0.6035533905932737308575
relative error = 5.194055e-02

division number = 4
integral value = 0.6284174365157310093721
relative error = 1.288420e-02

division number = 8
integral value = 0.6345731492255537453673
relative error = 3.214828e-03

division number = 16
integral value = 0.6361083632808495202937
relative error = 8.033195e-04

division number = 32
integral value = 0.6364919355013016755152
relative error = 2.008057e-04

division number = 64
integral value = 0.6365878141136421808355
relative error = 5.019991e-05

division number = 128
integral value = 0.6366117828642577913811
relative error = 1.254988e-05

division number = 256
integral value = 0.6366177749955105324986
relative error = 3.137465e-06

division number = 512
integral value = 0.6366192730247987041636
relative error = 7.843658e-07

division number = 1024
integral value = 0.6366196475319001457649
relative error = 1.960914e-07

division number = 2048
integral value = 0.6366197411586624888002
relative error = 4.902286e-08

division number = 4096
integral value = 0.6366197645653517422915
relative error = 1.225571e-08

division number = 8192
integral value = 0.6366197704170245552646
relative error = 3.063927e-09

division number = 16384
integral value = 0.6366197718799454507987
relative error = 7.659767e-10

division number = 32768
integral value = 0.6366197722456710117456
relative error = 1.914964e-10

division number = 65536
integral value = 0.6366197723371026517825
relative error = 4.787588e-11

division number = 131072
integral value = 0.6366197723599588131904
relative error = 1.197350e-11

division number = 262144
integral value = 0.6366197723656804585701
relative error = 2.985964e-12

シンプソン則

Simpson.exe
division number = 2
integral value = 0.6380711874576983078100
relative error = 2.279877e-03

division number = 4
integral value = 0.6367054518232166948621
relative error = 1.345850e-04

division number = 8
integral value = 0.6366250534621613610398
relative error = 8.295524e-06

division number = 16
integral value = 0.6366201012992814822766
relative error = 5.166847e-07

division number = 32
integral value = 0.6366197929081188755518
relative error = 3.226500e-08

division number = 64
integral value = 0.6366197736510885718531
relative error = 2.016128e-09

division number = 128
integral value = 0.6366197724477961061851
relative error = 1.260010e-10

division number = 256
integral value = 0.6366197723725945945006
relative error = 7.874735e-12

division number = 512
integral value = 0.6366197723678950204373
relative error = 4.926614e-13

division number = 1024
integral value = 0.6366197723676011444027
relative error = 3.104203e-14

division number = 2048
integral value = 0.6366197723675823816336
relative error = 1.569541e-15

division number = 4096
integral value = 0.6366197723675822706113
relative error = 1.395147e-15

division number = 8192
integral value = 0.6366197723675822706113
relative error = 1.395147e-15

division number = 16384
integral value = 0.6366197723675800501653
relative error = 2.092721e-15

division number = 32768
integral value = 0.6366197723675814934552
relative error = 1.743934e-16

division number = 65536
integral value = 0.6366197723675824926559
relative error = 1.743934e-15

division number = 131072
integral value = 0.6366197723675837139012
relative error = 3.662262e-15

division number = 262144
integral value = 0.6366197723675749431393
relative error = 1.011482e-14

ロンバーグ則

romberg.exe
0.6035533905932737308575
0.6284174365157310093721 0.6367054518232168058844
0.6345731492255537453673 0.6366250534621613610398 0.6366237772977001441177
0.6361083632808495202937 0.6366201012992814822766 0.6366200226935214701030 0.6366200079695835123417
0.6364919355013016755152 0.6366197929081190975964 0.6366197880130213260230 0.6366197870927056134960 0.6366197868767947687374
0.6365878141136421808355 0.6366197736510890159423 0.6366197733454218576910 0.6366197732879018689189 0.6366197732744074411215 0.6366197732710856538318
0.6366117828642577913811 0.6366197724477963282297 0.6366197724286963843809 0.6366197724251013712049 0.6366197724242580457954 0.6366197724240504340898 0.6366197724239986976968
0.6366177749955105324986 0.6366197723725948165452 0.6366197723714012157714 0.6366197723711765066312 0.6366197723711237710376 0.6366197723711107814282 0.6366197723711075617814 0.6366197723711066736030
0.6366192730247987041636 0.6366197723678946873704 0.6366197723678200803832 0.6366197723678059805508 0.6366197723678026498817 0.6366197723678018727256 0.6366197723678016506810 0.6366197723678016506810 0.6366197723678015396587
0.6366196475319001457649 0.6366197723676005892912 0.6366197723675959263545 0.6366197723675950381761 0.6366197723675948161315 0.6366197723675948161315 0.6366197723675948161315 0.6366197723675948161315 0.6366197723675948161315 0.6366197723675948161315
0.6366197411586624888002 0.6366197723675832698120 0.6366197723675829367451 0.6366197723675829367451 0.6366197723675829367451 0.6366197723675829367451 0.6366197723675829367451 0.6366197723675829367451 0.6366197723675829367451 0.6366197723675829367451 0.6366197723675829367451
0.6366197645653517422915 0.6366197723675814934552 0.6366197723675814934552 0.6366197723675814934552 0.6366197723675814934552 0.6366197723675814934552 0.6366197723675814934552 0.6366197723675814934552 0.6366197723675814934552 0.6366197723675814934552 0.6366197723675814934552 0.6366197723675814934552
0.6366197704170245552646 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890 0.6366197723675821595890
0.6366197718799454507987 0.6366197723675858233250 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473 0.6366197723675859343473
0.6366197722456710117456 0.6366197723675796060760 0.6366197723675794950537 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314 0.6366197723675793840314
0.6366197723371026517825 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207 0.6366197723675798281207
0.6366197723599588131904 0.6366197723675774966523 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300 0.6366197723675773856300
0.6366197723656804585701 0.6366197723675877107041 0.6366197723675878217264 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487 0.6366197723675879327487


division number = 1 integral value = 0.6035533905932737308575
relative error = 5.194055e-02
division number = 2 integral value = 0.6367054518232168058844
relative error = 1.345850e-04
division number = 3 integral value = 0.6366237772977001441177
relative error = 6.290930e-06
division number = 4 integral value = 0.6366200079695835123417
relative error = 3.700828e-07
division number = 5 integral value = 0.6366197868767947687374
relative error = 2.279102e-08
division number = 6 integral value = 0.6366197732710856538318
relative error = 1.419221e-09
division number = 7 integral value = 0.6366197724239986976968
relative error = 8.862011e-11
division number = 8 integral value = 0.6366197723711066736030
relative error = 5.537514e-12
division number = 9 integral value = 0.6366197723678015396587
relative error = 3.458222e-13
division number = 10 integral value = 0.6366197723675948161315
relative error = 2.110160e-14
division number = 11 integral value = 0.6366197723675829367451
relative error = 2.441508e-15
division number = 12 integral value = 0.6366197723675814934552
relative error = 1.743934e-16
division number = 13 integral value = 0.6366197723675821595890
relative error = 1.220754e-15
division number = 14 integral value = 0.6366197723675859343473
relative error = 7.150130e-15
division number = 15 integral value = 0.6366197723675793840314
relative error = 3.139082e-15
division number = 16 integral value = 0.6366197723675798281207
relative error = 2.441508e-15
division number = 17 integral value = 0.6366197723675773856300
relative error = 6.278163e-15
division number = 18 integral value = 0.6366197723675879327487
relative error = 1.028921e-14

結論

結論として、誤差が一番少ないものはシンプソン則とロンバーグ則であるが、ロンバーグ則は収束が早いが計算量が多く、シンプソン則は収束は遅いが計算量が少なくなっていることがわかった。

2
1
2

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
2
1