OdePack で常微分方程式を解く

インストール

準備として fpm をあらかじめインストールしておく (参考)。

git clone https://github.com/jacobwilliams/odepack.git
cd odepack/src/
make clean
make gfortran
make run
make test

自分の環境だと fpm build ではコンパイルできなかったので、Makefileを使った。
(external *** のように外部宣言している部分をコメントアウトすればうまくいくかも)
コンパイルすると M_odepack.o と m_odepack.mod が生成されるので、それぞれ /usr/local/lib、/usr/local/include などに移動しておく(任意).

計算準備

次の微分方程式を解く: $$ y'' -y' -2y = 18xe^{2x} $$

初期条件は $$ y(0) = \pi+\sqrt{2} = 4.5558062159628880 $$ $$ y'(0) = -2 + 2\pi - \sqrt{2} = 2.8689717448064913$$

理論解は以下の通り: $$ y' = (6x-2)e^{2x} + 2(3x^2-2x+\pi)e^{2x} - \sqrt{2}e^{-x} $$ $$ y = (3x^2-2x+\pi)e^{2x} + \sqrt{2}e^{-x} $$

ソースコード・コンパイル

lsoda を使って計算する。

program test
use m_odepack
implicit none

! odepack parameters
integer, parameter :: neq = 2
real(8) :: atol(neq)=1d-9, rtol=1d-9
real(8) :: y(neq), t, tout, rwork(200)
integer :: iopt=0
integer :: iout
integer :: istate=1, itask=1
integer :: itol=2
integer :: jt = 2
integer :: liw=100, lrw=200, iwork(100)

! function setting
integer, parameter :: N = 50
real(8), parameter :: tmin=0.d0, tmax=3.d0
real(8), parameter :: dt=(tmax-tmin)/dble(N-1)
real(8) :: time(N)
integer :: i

t = tmin                    ! initial time
y(1) = 4.5558062159628880d0 ! y(0)
y(2) = 2.8689717448064913d0 ! y'(0)

do i = 1, N
    if(i==1) time(i) = tmin
    if(i/=1) time(i) = time(i-1) + dt
end do

do i = 1, N
    tout = time(i)
    call dlsoda(fex,[neq],y,t,tout,itol,[rtol],atol,itask,istate,iopt,rwork,lrw,iwork,liw,jdum,jt)
    write (55,'(7ES12.4)') t, y(1), ftrue(t), y(1)-ftrue(t), y(2), dftrue(t), y(2)-dftrue(t)
    if ( istate<0 ) stop "istate < 0"
enddo

contains
! ================================================== !
subroutine fex( Neq, T, Y, Ydot ) 
integer, intent(in) :: neq
real(8), intent(in) :: t, y(neq)
real(8), intent(out) :: ydot(neq)
   Ydot(1) = Y(2)
   Ydot(2) = Y(2) + 2.d0*Y(1) + 18.d0*t*exp(2.d0*T)
end subroutine fex
! ================================================== !
subroutine jdum()
implicit none
end subroutine jdum
! ================================================== !
real(8) function ftrue(x)

real(8), intent(in) :: x
real(8) :: k1 = sqrt(2.d0), k2 = acos(-1.d0)

ftrue = (3.d0*x**2 -2.d0*x + k2)*exp(2.d0*x) + k1*exp(-x)

end function ftrue
! ================================================== !
real(8) function dftrue(x)
    real(8), intent(in) :: x
    real(8) :: k1, k2
    k1 = sqrt(2.d0)
    k2 = acos(-1.d0)
    dftrue = (6.d0*x - 2.d0)*exp(2.d0*x) + 2.d0*(3.d0*x**2 - 2.d0*x + k2)*exp(2.d0*x) - k1*exp(-x)
end function dftrue
! ================================================== !
end program test

• neq は微分方程式の階数 (変数の数)

• atol は絶対誤差 (abs)、rtolは相対誤差 (rel)

• その他は odepack の手続き変数

• 解く時間は自分で設定した. 最初に t = tmin として、ループでは t -->> tout になるように時間が更新される。

• $y(1)$ を 求める $y$、$y(2)$ を $y'$ とすると、$ydot(1)=y' = y(2)$、$ydot(2)=y(2)+2y(1)+18xe^{2x}$ と書けるので、これを解く。dotは微分を意味。

• $x<0$ も解けるが、初期値が必要. 今は理論解がわかっているのですぐに試せる.

コンパイル

gfortran /usr/local/lib/M_odepack.o test.f90 #-I/usr/local/include/

実行結果

$ head fort.55
  0.0000E+00  4.5558E+00  4.5558E+00  0.0000E+00  2.8690E+00  2.8690E+00  0.0000E+00
  6.1224E-02  4.7554E+00  4.7554E+00  1.8396E-09  3.6747E+00  3.6747E+00  5.7735E-09
  1.2245E-01  5.0092E+00  5.0092E+00  2.3605E-09  4.6483E+00  4.6483E+00  6.6694E-09
  1.8367E-01  5.3288E+00  5.3288E+00  3.5579E-09  5.8302E+00  5.8302E+00  9.3923E-09
  2.4490E-01  5.7283E+00  5.7283E+00  6.7996E-09  7.2697E+00  7.2697E+00  1.7355E-08
  3.0612E-01  6.2254E+00  6.2254E+00  6.3543E-09  9.0258E+00  9.0258E+00  1.5290E-08
  3.6735E-01  6.8415E+00  6.8415E+00  7.4869E-09  1.1170E+01  1.1170E+01  1.7392E-08
  4.2857E-01  7.6028E+00  7.6028E+00  9.2457E-09  1.3788E+01  1.3788E+01  2.1080E-08
  4.8980E-01  8.5416E+00  8.5416E+00  1.1599E-08  1.6984E+01  1.6984E+01  2.6115E-08
  5.5102E-01  9.6968E+00  9.6968E+00  1.6157E-08  2.0880E+01  2.0880E+01  3.6313E-08
  時刻 t      計算した y   理論解 yt    誤差 y-yt   計算したy'   理論解 y't  誤差 y'-y't

誤差は自分で設定した tol = 1d-9 のオーダーになっている。
例えば変数を減らす (ex: y' = ky) ときは neq = 1、増やすときは neq = 3 のようにする

参考

chatgptさんによると、

lsoda という解法について教えて！


ChatGPT:
いい質問です！LSODA は 常微分方程式 (ODE) を数値的に解くための非常に有名で賢いアルゴリズムのひとつです。
Fortran で最初に実装されましたが、今では Python（SciPy）、MATLAB、Julia など色々な言語で使えるようになっています。

🔧 LSODA の特徴
LSODA は 自動的に「剛性のある方程式かどうか」を判断して、適切な手法を切り替えてくれるハイブリッドなソルバ です。

剛性でない場合（stiffでない）：
→ Adams method（マルチステップ法）を使用

剛性がある場合（stiff）：
→ BDF法（後退差分公式） を使用

これによって、ユーザーが「剛性かどうか」を気にせずに使える、便利なソルバになっています。

📌 LSODA の使いどころ
方程式が剛性かどうか分からないとき

硬い部分とそうでない部分が混在しているような現象（例：化学反応、プラズマ、流体など）

長時間の安定したシミュレーションが必要なとき

🔍 Python（SciPy）での使い方例
python
コピーする
編集する
from scipy.integrate import solve_ivp
import numpy as np
import matplotlib.pyplot as plt

# 例：dy/dt = -y, 初期値 y(0) = 1
def dydt(t, y):
    return -y

sol = solve_ivp(dydt, [0, 10], [1], method='LSODA')

plt.plot(sol.t, sol.y[0])
plt.xlabel("t")
plt.ylabel("y")
plt.title("dy/dt = -y を LSODA で解いた結果")
plt.grid()
plt.show()
🧠 補足：剛性（Stiffness）って？
剛性のある方程式とは、急激に変化する成分と緩やかに変化する成分が混在していて、普通の数値解法だとすごく小さなステップが必要になってしまうような ODE のこと。

LSODA のような「剛性対応ソルバ」がないと、解くのに時間がかかりすぎたり、不安定になったりします。

もし Fortran で使いたい場合や、実際に解きたい方程式があるなら、それに合わせてコード例を出すこともできるよ！どうする？

だそうです。。