lsoda の使い方

R で lsoda を使って微分方程式を数値的に解く。

R には lsoda というパッケージがある。これを導入すると、常微分方程式の数値積分を行うことができる。
導入の仕方などについては、「Windows に Eclipse + R を導入する」を参照のこと。

使い方をいろいろ書くよりも、サンプルを示した方が簡単だと思うので、いくつかのサンプルを示す。

• 自由落下の例

         
              library(odesolve)
              
              #
              #  時間変化を記述する関数(微分方程式の右辺を与える関数)
              #
              #      引数   : t(現在時刻), y(現在の変数), parms(パラメタ)
              #      返り値 : y の時間微分を list(c( )) の中に並べる。
              #
              #      今回の場合、
              #        y[1] : 座標
              #        y[2] : 速度
              #
              freefall <- function(t, y, parms){
                y1dot <- y[2]
                y2dot <- - parms['g']
                
                return( list(c(y1dot, y2dot)) )
              }
              
              #
              #  パラメタ
              #      上記関数の中で参照するために、名前をつけておく。
              #
              parms <- c( g=9.8, m=1 )
              
              #
              #  時刻
              #      計算結果を表示したい時刻の値を配列で与える。
              #
              times <- seq(0,10,0.01)
              
              #
              #  初期値
              #         ここで変数に名前をつけておくこともできる。
              #
              y     <- c( y=0, vy=0 )
              
              #
              #  計算と結果の記録
              #
              result <- lsoda(y, times, freefall, parms)
              
              #
              #  表示
              #
              plot(0, type='n',
                   xlim=c(min(times), max(times) ),
                   ylim=c(min(result[,'y']), max(result[,'y'])))
              
              # #
              # #  本来ならば
              # points( result[,1], result[,2], pch='.')
              # #  などで表示がすむところを、あえて時間をかけて表示させる。
              # #  ちょっとアニメーション気分で
              
              for( i in 1:length(result[,2]) ){
                points( result[i,'time'], result[i,'y'], pch='.' , col='red')
              }
         

• 地球 - 月　系

         
              library(odesolve)
              
              #
              #  時間変化を記述する関数(微分方程式の右辺を与える関数)
              #
              #      引数   : t(現在時刻), y(現在の変数), parms(パラメタ)
              #      返り値 : y の時間微分を list(c( )) の中に並べる。
              #
              #      今回の場合、
              #        y[1] : 地球の x 座標
              #        y[2] : 地球の y 座標
              #        y[3] : 地球の 速度 x成分
              #        y[4] : 地球の 速度 y成分
              #
              #        y[5] : 月の x 座標
              #        y[6] : 月の y 座標
              #        y[7] : 月の 速度 x成分
              #        y[8] : 月の 速度 y成分
              #
              Earth <- function(t, y, parms){
                const     <- parms['G'] * parms['massE'] * parms['massM']
                distance2 <- (y[1] - y[5])^2 + (y[2] - y[6])^2
                distance  <- sqrt(distance2)
                
                y1dot <- y[3]
                y2dot <- y[4]
                y3dot <- const / distance2 * ( y[5] - y[1] ) / distance / parms['massE']
                y4dot <- const / distance2 * ( y[6] - y[2] ) / distance / parms['massE']
              
                y5dot <- y[7]
                y6dot <- y[8]
                y7dot <- const / distance2 * ( y[1] - y[5] ) / distance / parms['massM']
                y8dot <- const / distance2 * ( y[2] - y[6] ) / distance / parms['massM']
              
                return( list(c(y1dot, y2dot, y3dot, y4dot,
                               y5dot, y6dot, y7dot, y8dot)) )
              }
              
              #
              #  パラメタ
              #      上記関数の中で参照するために、名前をつけておく。
              #
              parms <- c( G=6.67300E-11, massE=5.9742E24, massM=7.347673E22 )
              
              #
              #  時刻
              #      計算結果を表示したい時刻の値を配列で与える。
              #      1 時間ごとに 90 日分
              #
              times <- seq(0, 60*60*24*90, 60*60)
              
              #
              #  初期値
              #         ここで変数に名前をつけておくこともできる。
              #
              #         地球 x   y     vx     vy          月 x   y     vx     vy
              y     <- c( Ex=0, Ey=0, Evx=0, Evy=0, Mx=3.844E8, My=0, Mvx=0, Mvy=1.023144e3)
              
              #
              #  計算と結果の記録
              #
              result <- lsoda(y, times, Earth, parms)
              
              #
              #  表示
              #
              plot(0, type='n',
                   xlim=c(-max(result[,'My']), max(result[,'My'])),
                   ylim=c(-max(result[,'My']), max(result[,'My'])))
              
              # #
              # #  本来ならば、
              # points( result[,'Ex'], result[,'Ey'], pch='.', col='blue')
              # points( result[,'Mx'], result[,'My'], pch='.', col='orange')
              # #  で済むところを、ちょっとアニメーション気分で。
              # #
              for( i in 1:length(result[,2]) ){
                points( result[i,'Ex'], result[i,'Ey'], pch='.' , col='blue')
                points( result[i,'Mx'], result[i,'My'], pch='.' , col='orange')
              }