微分方程式の数値解析 #10

有名どころでLorenz方程式を解いてみる。
方程式のパラメータはLorenzが最初の論文で示したものを使用する。

program ode
  use runge_kutta_4
  implicit none
  double precision, parameter :: a = 10, b = 28, c = 8d0/3    ! Lorenz's parameters
  double precision, parameter :: t0 = 0                       ! initial time
  double precision, parameter :: x0(3) = (/1d-1, 1d-1, 1d-1/) ! initial coordinate at t0
  double precision :: t
  double precision :: x(size(x0))
  double precision, parameter :: udiv = 1         ! division unit
  integer, parameter :: ndiv = 1000               ! number of divisions per udiv
  double precision, parameter :: h = udiv / ndiv  ! delta-t
  integer, parameter :: niter = 50 / udiv * ndiv  ! number of iterations
  integer, parameter :: pdiv = ndiv / 100         ! print interval
  integer :: n
  type (func_type), dimension(size(x0)) :: fns
  fns(1)%f => f1
  fns(2)%f => f2
  fns(3)%f => f3
  call set_func(fns)
  x = x0
  do n = 0, niter
    t = t0 + n * h
    if (mod(n, pdiv) == 0) print '(4F13.8)', t, x(1), x(2), x(3)
    x = next(t, x, h)
  end do
contains
  double precision function f1(t, x)
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f1 = a * (x(2) - x(1))
  end function f1
  double precision function f2(t, x)
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f2 = x(1) * (b - x(3)) - x(2)
  end function f2
  double precision function f3(t, x)
    double precision, intent(in) :: t
    double precision, dimension(:), intent(in) :: x
    f3 = x(1) * x(2) - c * x(3)
  end function f3
end program ode

結果をgnuplotで描画する。

set xlabel "X"
set ylabel "Y"
set zlabel "Z"
set terminal png
set output "lorenz.png"
splot 'lorenz.txt' using 2:3:4 with lines
set output