Julia(1.4.1)~常微分方程式の数値積分~
やること
Julia言語で数値積分を行う．
数値積分に使うものとしてDifferensialEquations.jlとSundials.jlがある模様だが，DifferentialEquationsを使用する． 基本的には以下の公式サイトに乗っている．
Home · DifferentialEquations.jl
環境
- JupyterLab
- "julia" => v"1.4.1"
- "DifferentialEquations" => v"6.13.0"
- "PyPlot" => v"2.9.0"
常微分方程式（一変数，式が一つの場合）
次の数式を積分する．
df/dt = cos(t)
積分したら，f = sin(t)となるはず．
モジュールのインストール
DifferentialEquationsがない場合は以下のコードをJupyterLabのセルで実行する．
import Pkg Pkg.update() Pkg.add("DifferentialEquations")
モジュールの定義とプリコンパイル
using DifferentialEquations #結果表示用 using PyPlot
数値積分のコード
まずは積分の準備を行う
#関数の定義 f(u,p,t) = cos(t) #初期値の定義 u0 = 0.0 #積分区間の指定 tspan = (0.0,4*pi) #積分問題の定義 prob = ODEProblem(f,u0,tspan)
ポイント
f(u,p,t)について，uは変数，tは時間，pは外部からの入力を与えるもの．- 初期値
u0は整数にすると計算が回らない．~.0のようにFloat型になるようにする． - 積分区間の指定も同様．始まりと終わりを指定するが，両者Floatになるようにする．
次に，積分を実行する． ただ単に積分するだけなら以下のコード
#セミコロンをつけると結果をprintしない
sol = solve(prob);
結果のプロット
# pygui(true)とすることでインタラクティブなwindow出力 # pygui(false) #sol.tが時間，sol.uに出力が格納される plot(sol.t,sol.u) #画像の保存 savefig("sin.png")
以下のプロットが出力される．
流石に画像が荒いので，出力される間隔を0.1ごとにする．
積分の実行コードに saveat=0.1 を追加．
sol = solve(prob,saveat=0.1);
結果をプロットすると以下のようになめらかな出力が得られる．
solverの変更
solverを指定することができる．たとえば，Tsit5を指定するときは以下．
なにも指定しないときはTsit5が適用される．
sol = solve(prob,Tsit5(),saveat=0.1);
その他のソルバーについては以下参照
ODE Solvers · DifferentialEquations.jl
matlabとの対応は以下．
solverの対応と種類
| solver | 対応するMATLABのsolver | 説明 |
|---|---|---|
| BS3() | ode23 | 低精度，高速 |
| DP5() | ode45 | 中精度 |
| Tsit5() | ode45 | 中精度 |
| VCABM() | ode113 | 低～高精度 |
| Vern9() | - | 高精度 |
ソルバーによる精度と実行時間の検証はまた今度行いたい．
相対誤差と絶対誤差の指定
相対誤差と絶対誤差を指定することで，より精度の高い出力を得られる．誤差をどこまで小さく指定できるかはsolverによる（Tsit5で1e-8程度）
相対誤差はreltol =，絶対誤差はabstol =で指定する．
sol = solve(prob,Tsit5(),reltol = 1e-8,abstol = 1e-8,saveat=0.1);
常微分方程式（変数２つの場合）
2つ以上の変数を有する方程式は，関数宣言の部分が少し異なる．具体的には関数名に!をつける．
これはjuliaの関数で，配列を引数に取る場合!をつける仕様によるもの．
例
以下の式を積分する
dx/dt = x-y dy/dt = y
コード
#引数にduを追加する．u[1] = x,u[2]=yとなる #function名f2に!をつける function f2!(du,u,p,t) du[1] = u[1]-u[2] du[2] = u[2] end #初期値，u0 = [x0,y0] u0 = [1.0,1.0] tspan = (0.0,4*pi) #ここでも関数名に!をつけるのを忘れない prob = ODEProblem(f2!,u0,tspan) sol2 = solve(prob)
以上
