2.7 計算の精度と安定性：問題 2.3 オイラー法の計算精度と安定性 a, b, c

今回は計算精度をマジメに考えようって回です。　対象とする物理系はニュートンの冷却の法則にしたがう系です。　この物理系については以前の記事参照。　今回は問題2.3のうち、a, b, c をやっていきます。　d, e は次回に。

問題

a. ニュートンの冷却の法則の解析解が
　　 $T(t) = T_s - (T_s - T_0) e^{-rt}$
となることを示せ。
b. プログラム cool を使って、Δt = 0.1, 0.05, 0.025, 0.01, 0.005 について 1min の温度を計算せよ。　解析解と、得られた数値解の差を、Δt の関数として示す表を作れ。　オイラー法の次数はいくらか。
c. t = 1 で 0.1% の（3桁の）精度を得るには Δt の値をどうとったらよいか。　t = 4 で 0.1% の（3桁の）精度を得るには Δt の値をどうとったらよいか。

bに対して: t の値を一定にしたときに解析解と数値解の差が (Δt)ⁿ に比例するなら、その数値計算法は n 次であるという。
c に対して: 数値計算の精度を決定する1つの方法は、より小さい分割幅にして計算を繰り返し実行し、それらの結果を比較することである。　2つの計算が p 桁まで一致するなら、計算結果は p 桁まで正しいと見なしてよいだろう。

a.

ニュートンの冷却の法則を表す微分方程式

　　 $\frac{dT}{dt} = -r(T - T_s)$

は最も簡単な微分方程式の部類なので、解くのは簡単かと。　一応こちらで導いてみました。

b.

次はいろいろな時間分割幅 dt でシミュレーションを行います。　並行シミュレーションでいきましょう：

@GrabResolver('https://github.com/waman/groops-core/tree/master/repo')
@Grab('org.waman.groops:groops-core:0.1-beta')

import org.waman.groops.builder.SimulationBuilder
import org.waman.groops.simulation.system.*

double T0 = 82.3d, T_room = 17.0d, r = 0.03d
def tmax = 1.0

// シミュレーションする dt のリスト
def dtList = [0.1, 0.05, 0.025, 0.01, 0.005]

// 解析解
def analyticSolution =  { BigDecimal t -> T_room - (T_room - T0)*Math.exp(-r*t.doubleValue()) }
double analyticT = analyticSolution(tmax)

new SimulationBuilder().parallel{

    def pngout = png(fileName:'problem2_3b.png', overwrite:true, outputTiming:FINAL){
        chart(title:'problem 2.3 b'){
            point(label:'Simulation', x:'dt', y:{ it.T - analyticT })
        }
    }

    for(dt in dtList){
        simulation{
            system(new CoolSystem(), dt:dt, T_room:T_room, r:r)
            init(time:0.0, T:T0)
            parameterObservableSet()

            continueWhile{ it.time <= tmax }
            dataOutputter(pngout['Simulation'])
        }
    }
}.simulate()

class CoolSystem extends AnnotationObservableSystem{

    @State BigDecimal time
    @State double T

    @Parameter BigDecimal dt
    @Parameter double T_room
    @Parameter double r
    
    @Override
    void evolve(Map params){
        T += -r * (T - T_room)*dt.doubleValue()
        time += dt
    }
}

problem2_3b.groovy

結果は下図のようになります：

結構キレイに直線になりますね。　回帰直線とか出してませんけど。　キチンと解析するなら、データをファイル等に書き出して別途解析を行った方がいいでしょうね。　まぁ、ともかくオイラー法の次数は1次だと分かりました。

c.

次もいろいろな dt の値のシミュレーションを並行して行いましょう。　ただし、今回の精度判定は、シミュレーション結果を書き出して p 桁の値が一致しているかどうかを自分で判定することにします。　下記のコードは一部：

def dtList = [0.1, 0.05, 0.025, 0.01, 0.005, 0.0025, 0.001]

double T0 = 82.3d, T_room = 17.0d
double r = 0.03d
def tmax = 1.0 // or 4.0

new SimulationBuilder().parallel{
    for(dt in dtList){
        simulation{
            system(new CoolSystem(), dt:dt, T_room:T_room, r:r)
            init(time:0.0, T:T0)
            parameterObservableSet()

            continueWhile{ it.time <= tmax }
            console(dataEntries:['dt', 'T'], outputTiming:FINAL)
        }
    }
}.simulate()

problem2_3c.groovy

結果を表にすると以下のようになります:

dt	t = 1 での温度	t = 4 での温度
0.1	80.17713432665997	74.73174332091298
0.05	80.2736131031035	74.82382107228912
0.025	80.32185303946834	74.86986208399091
0.01	80.3507971764579	74.89748737364134
0.005	80.36044525129657	74.90669591722178
0.0025	80.36526929418437	74.91130021032322
0.001	80.36816372166693	74.91406279300254

3桁の精度になるのは

t=1 の場合： dt = 0.01
t=4 の場合： dt = 0.025

でした。　ついでに4桁の精度も調べると

t=1 の場合： dt = 0.025
t=4 の場合： dt = 0.001

です。

3桁の精度の場合、時間が多く経過した t=4 の場合の方が大きい dt で精度が達成されています。　これは精度を桁数の数字だけで評価しているせいでしょう。　この精度の評価方法では 0.09 に対して 0.1 より 0.01 の方が近いとされてしまうので、ちょっと直感に反した評価結果になることがあります。

他にも、今回の場合は dt の細分化方法（0.1 → 0.05 → 0.025 ・・・）が一様ではないというのも注意が必要です。　きちんとやるには、常に分割数を2倍にするなどの方がよいでしょう。