Mercurial > op > rk_sakamoti
changeset 0:5bd0f2a15d2a
import from https://qiita.com/sakamoti/items/de851e0d07aeef5be310
| author | "uncorrelated zombie" <uncorrelated@yahoo.co.jp> |
|---|---|
| date | Tue, 23 Jul 2024 02:30:57 +0900 |
| parents | |
| children | 954faa9f6837 |
| files | .hgignore rkmain.f90 rungekutta.f90 |
| diffstat | 3 files changed, 135 insertions(+), 0 deletions(-) [+] |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/.hgignore Tue Jul 23 02:30:57 2024 +0900 @@ -0,0 +1,4 @@ +syntax: glob +*.mod +*.o +
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rkmain.f90 Tue Jul 23 02:30:57 2024 +0900 @@ -0,0 +1,31 @@ +program rkmain + use rungekutta !ルンゲクッタ解法をまとめたモジュール + implicit none + real(kind=dpkind),dimension(3) :: x !状態変数 + type(rk) :: rktype !ルンゲクッタ法ソルバに関連する型 + + !微分方程式の右辺を計算する関数、解法、時間刻みの指定 + call rktype%setup(rhside,"rk4",1d-5) + !初期値の指定 + x=[0.5e0,0.5e0,0.5e0] + !時刻0?50秒まで、0.01秒毎に標準出力へ書き出し。 + !ただし、内部的には刻み幅1e-5秒としている。 + do + if(rktype%time>5e1)exit + call rktype%solver(rktype%time+1e-2,x) !積分実行 + print*,rktype%time,x + enddo + !======================================== + contains + !微分方程式の右辺を定義 + function rhside(self,time,wk) + class(rk),intent(inout) :: self !parameter + real(kind=dpkind) :: time + real(kind=dpkind),dimension(:) :: wk + real(kind=dpkind),dimension(size(wk)) :: rhside + real(kind=dpkind),parameter :: a=10e0,b=28e0,c=8e0/3e0 + rhside(1)=a*(wk(2)-wk(1)) + rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3) + rhside(3)=wk(1)*wk(2)-c*wk(3) + end function +end program
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/rungekutta.f90 Tue Jul 23 02:30:57 2024 +0900 @@ -0,0 +1,100 @@ +module rungekutta + use,intrinsic :: iso_fortran_env + private + !--- + integer,parameter,public :: dpkind=real64 + real(kind=dpkind),public,parameter :: PI=3.1415926535897932384626433_dpkind + !--- + type,public :: rk + real(kind=dpkind) :: time=0e0 !<險育ｮ鈴幕蟋九°繧峨ｮ邱冗ｵ碁℃譎る俣 + real(kind=dpkind) :: dt=1e-4 !<譎る俣蛻ｻ縺ｿ + procedure(rhside),pass,pointer :: fun => null() + !<蠕ｮ蛻譁ｹ遞句ｼ上ｮ蜿ｳ霎ｺ繧定ｨ育ｮ励☆繧矩未謨ｰ + procedure(rksol),pass,pointer :: solver => null() + !<蠕ｮ蛻譁ｹ遞句ｼ上ｮ隗｣豕 + contains + final :: final_rk !<讒矩菴薙ｮ繝輔ぃ繧､繝翫Λ繧､繧ｺ蜃ｦ逅 + procedure :: setup !<謨ｰ蛟､遨榊縺ｮ蛻晄悄蛟､縲∬ｧ｣豕輔ｒ繧ｻ繝繝医☆繧区焔邯壹″ + end type + !------------- + interface + !--- + ! Right hand side function interface + function rhside(self,time,wk) + import :: dpkind,rk + class(rk),intent(inout) :: self !parameter + real(kind=dpkind) :: time + real(kind=dpkind),dimension(:) :: wk + real(kind=dpkind),dimension(size(wk)) :: rhside + end function + !--- + ! Solver interface + subroutine rksol(self,tnext,wk) + import :: dpkind,rk + class(rk),intent(inout) :: self + real(kind=dpkind),intent(in) :: tnext !output time + real(kind=dpkind),dimension(:) :: wk + end subroutine + end interface + !============================-- + contains + !--- + pure elemental subroutine final_rk(self) + type(rk),intent(inout) :: self + if(associated(self%fun))nullify(self%fun) + if(associated(self%solver))nullify(self%solver) + end subroutine + !--- + subroutine setup(self,fun,solver,dt) + class(rk),intent(inout) :: self + procedure(rhside) :: fun !蠕ｮ蛻譁ｹ遞句ｼ上ｮ蜿ｳ霎ｺ繧定ｨ育ｮ励☆繧矩未謨ｰ + character(len=*),intent(in) :: solver !繧ｽ繝ｫ繝舌ｼ縺ｮ驕ｸ謚 + real(kind=dpkind),intent(in),optional :: dt !譎る俣蛻ｻ縺ｿ蟷 + !譎る俣蛻ｻ縺ｿ縺ｮ繧ｻ繝繝 + if(present(dt))then + self%dt=dt + endif + self%fun=>fun + !繧ｽ繝ｫ繝舌ｼ縺ｮ驕ｸ謚 + select case(solver) + case("rk4") + !繝ｫ繝ｳ繧ｲ繧ｯ繝繧ｿ豕 + self%solver=>rk4 + case("rkf45") + !繝ｫ繝ｳ繧ｲ繧ｯ繝繧ｿ繝ｻ繝輔ぉ繝ｫ繝吶Ν繧ｰ豕(蛻ｻ縺ｿ蟷閾ｪ蜍戊ｪｿ謨ｴ蝙) + !self%solver=>rkf45 + case default + write(error_unit,*) "solver must one of the 'rk4, rkf45'" + endselect + end subroutine + !--- + subroutine rk4(self,tn,wk) + class(rk),intent(inout) :: self + real(kind=dpkind),intent(in) :: tn !谺｡縺ｫ蜃ｺ蜉帙＠縺溘＞譎る俣 + real(kind=dpkind),dimension(:) :: wk + real(kind=dpkind),dimension(size(wk)) :: k1,k2,k3,k4 + real(kind=dpkind) :: dt_tmp + real(kind=dpkind) :: dthalf + associate(time=>self%time,dt=>self%dt) + dt_tmp=dt + do + if(time+dt>tn)dt_tmp=tn-time + dthalf=dt_tmp/2e0 + !!$OMP workshare + k1=self%fun(time,wk) + k2=self%fun(time+dthalf,wk+k1*dthalf) + k3=self%fun(time+dthalf,wk+k2*dthalf) + k4=self%fun(time+dt_tmp,wk+k3*dt_tmp) + !迥ｶ諷句､画焚繧呈峩譁ｰ + wk=wk+dt_tmp*(k1+2e0*k2+2e0*k3+k4)/6e0 + !!$OMP end workshare + !1繧ｹ繝繝繝玲ｬ｡縺ｮ隗｣縺瑚ｨ育ｮ励〒縺阪◆縺ｮ縺ｧ縲∵凾髢薙ｒ騾ｲ繧√ｋ + time=time+dt_tmp + !譎る俣縺ｮ險育ｮ礼ｯ蝗ｲ縺荊n繧定ｶ縺医◆繧芽ｨ育ｮ礼ｵゆｺ + if(time>=tn)then + exit + endif + enddo + end associate + end subroutine +end module