微分方程式の数値解析 #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