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

今回は、前回に続き計算精度の問題をやっていきます。　対象とする物理系はガラッと変わって RC 回路になります。　といっても、この系が従う法則もニュートンの冷却の法則と同じ形の微分方程式なので、結局やってることは変わりません。

問題

d. 微分方程式
　　 $R\frac{dQ}{dt} = V - \frac{Q}{C}$
を考える。　ただし ${ t = 0 }$ で ${ Q = 0 }$ とする。　この方程式は加えられた電圧が V の時の RC 回路のコンデンサーの充電過程を表す。　 ${ t }$ については秒を単位にとり、 ${ R = 2000\,[\Omega],\,C = 10^{-6}\,[\textrm{F}],\,V = 10\,[\textrm{V}] }$ として、この微分方程式をオイラー法で解くプログラムを書け。　t = 0.005s では3桁の精度を得るには Δt の値をどうとる必要があるか。
e. Δt = 0.005, 0.0025, 0.001 としたときの t = 0.005s における RC 回路の数値解の性質を調べよ。

RC 回路の解析解はこちらで求めてます。　結果は（ ${ t = 0 }$ のとき ${ Q = 0 }$ の初期条件のもとで）

　　 $Q(t) = CV\left(1- e^{-\frac{t}{RC}}\right)$

となります。

d.

では RC 回路のシミュレーション。　と言っても、やってることは問題 2.3 c と同じです。　物理系のドメイン分析だけやっておきましょう：

状態変数
- 経過時間 ${ t }$ [s] ： t (BigDecimal)
- コンデンサーに充電されている電荷 ${ Q }$ [C] ： q (double)
パラメータ
- 時間間隔 ${ \Delta t }$ [s] ： dt (BigDecimal)
- 電圧 ${ V }$ [V] ： v (double)
- 抵抗 ${ R }$ ${ \Omega }$ ： r (double)
- 静電容量 ${ C }$ [F] ： c (double)

これを踏まえて、シミュレーション・スクリプトは以下のようになります：

@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.*

def dtList = [5.0E-4, 2.5E-4, 1.0E-4, 5.0E-5, 2.5E-5, 1.0E-5, 0.5E-5]
def tmax = 5.0E-3

new SimulationBuilder().parallel{

    for(dt in dtList){
        simulation{
            system(new RCCircuitSystem(), dt:dt, r:2000d, c:1E-6d, v:10d)
            init(t:0.0, q:0d)
            parameterObservableSet()

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

class RCCircuitSystem extends AnnotationObservableSystem{

    @State BigDecimal t
    @State double q

    @Parameter BigDecimal dt
    @Parameter double v, r, c
    
    @Override
    void evolve(Map params){
        double change = (v - q/c)/r
        q += change * dt.doubleValue()
        t += dt
    }
}

problem2_3d.groovy

シミュレーション結果は下表のようになります：

dt	t = 0.005s における結果
0.0005	9.577648639678956E-6
0.00025	9.39442336072811E-6
0.0001	9.269022734871224E-6
0.00005	9.22470642384714E-6
0.000025	9.20208254807407E-6
0.00001	9.188359978669231E-6
0.000005	9.183761147410228E-6

よって、3桁の精度を得るには dt = 0.000005 以下にとる必要があります。

e.

こちらは問題 2.3 b と同じような問題です。　上記のシミュレーションと同じような箇所（import 文、RCCircuitSystem クラスの定義）は省略して、シミュレーション・コードは以下のようになります：

double r = 2000d, c = 1E-6d, v = 10d
def tmax = 5.0E-3

def analyticSolution = { BigDecimal t -> c*v*(1d - Math.exp(-t.doubleValue()/(r*c)))}
double analyticQ = analyticSolution(tmax)

def dtList = [5.0E-4, 2.5E-4, 1.0E-4]

new SimulationBuilder().parallel{

    def pngout = png(fileName:'problem2_3e.png', overwrite:true, outputTiming:FINAL){
        chart(title:'problem 2.3 e', rangeLabel:'Q', range:[0.9E-5d, 1.0E-5d]){
            point(label:'Simulation', x:'dt', y:'q')
            marker(label:'Analytic', y:analyticQ, color:java.awt.Color.BLUE)
        }
    }

    for(dt in dtList){
        simulation{
            system(new RCCircuitSystem(), dt:dt, r:r, c:c, v:v)
            init(t:0.0, q:0d)
            parameterObservableSet()

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