Mercurial > op > rk_sakamoti

changeset 2:206e90e91191

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
make rkmod.f90
author "uncorrelated zombie" <uncorrelated@yahoo.co.jp>
date Sun, 21 Jul 2024 03:56:49 +0900
parents 954faa9f6837
children 2347cff808b0
files rkmain.f90 rkmod.f90
diffstat 2 files changed, 46 insertions(+), 43 deletions(-) [+]
line wrap: on
line diff
--- a/rkmain.f90	Sun Jul 21 02:08:58 2024 +0900
+++ b/rkmain.f90	Sun Jul 21 03:56:49 2024 +0900
@@ -1,49 +1,13 @@
 program rkmain
-  use rungekutta !ルンゲクッタ解法をまとめたモジュール
+  use rkmod
   implicit none
-  real(kind=dpkind),dimension(3) :: x !状態変数
-  type(rk) :: rktype  !ルンゲクッタ法ソルバに関連する型
-  real(kind=dpkind),allocatable,dimension(:,:) :: r
+  real(kind=dpkind),dimension(3) :: x
   real(kind=dpkind) :: max_t
-
-  !初期値の指定
-  x = [0.5e0,0.5e0,0.5e0]
-  max_t = 5e1
-  call calc(x, max_t, r)
-  print *, r;
+  real(kind=dpkind),allocatable,dimension(:,:) :: r
 
-  !========================================
-  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
-    !Rから呼び出すサブルーチン
-    subroutine calc(x, max_t, r)
-      implicit none
-      real(kind=dpkind),dimension(3),intent(in) :: x
-      real(kind=dpkind),intent(in) :: max_t
-      real(kind=dpkind),allocatable,dimension(:,:),intent(out) :: r
-      real(kind=dpkind) :: zero = 0
-      integer i, max_i
-      !微分方程式の右辺を計算する関数、解法、時間刻みの指定
-      call rktype%setup(rhside,"rk4",1d-5) 
-      max_i = int(max_t / 0.01)
-      allocate(r(0:max_i, 4))
-      !時刻max_t秒まで、0.01秒毎に配列に記録。
-      !ただし、内部的には刻み幅1e-5秒としている。
-      r(0, :) = [zero, x] ! 開始点
-      do i = 1, max_i
-        call rktype%solver(rktype%time+1e-2,x) !積分実行
-        r(i, :) = [rktype%time, x]
-      end do
-    end subroutine
+  x = [0.5d0,0.5d0,0.5d0]
+  max_t = 5d1
+  call calc(x, max_t, r)
+  print *, r
 end program
 
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rkmod.f90	Sun Jul 21 03:56:49 2024 +0900
@@ -0,0 +1,39 @@
+module rkmod
+  use rungekutta !繝ｫ繝ｳ繧ｲ繧ｯ繝繧ｿ隗｣豕輔ｒ縺ｾ縺ｨ繧√◆繝｢繧ｸ繝･繝ｼ繝ｫ
+  implicit none
+  type(rk) :: rktype  !繝ｫ繝ｳ繧ｲ繧ｯ繝繧ｿ豕輔た繝ｫ繝舌↓髢｢騾｣縺吶ｋ蝙
+  !========================================
+  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
+    !R縺九ｉ蜻ｼ縺ｳ蜃ｺ縺吶し繝悶Ν繝ｼ繝√Φ
+    subroutine calc(x, max_t, r)
+      implicit none
+      real(kind=dpkind),dimension(3),intent(in) :: x !迥ｶ諷句､画焚
+      real(kind=dpkind),intent(in) :: max_t
+      real(kind=dpkind),allocatable,dimension(:,:),intent(out) :: r
+      real(kind=dpkind) :: zero = 0
+      integer i, max_i
+      !蠕ｮ蛻譁ｹ遞句ｼ上ｮ蜿ｳ霎ｺ繧定ｨ育ｮ励☆繧矩未謨ｰ縲∬ｧ｣豕輔∵凾髢灘綾縺ｿ縺ｮ謖螳
+      call rktype%setup(rhside,"rk4",1d-5) 
+      max_i = int(max_t / 0.01)
+      allocate(r(0:max_i, 4))
+      !譎ょ綾max_t遘偵∪縺ｧ縲0.01遘呈ｯ弱↓驟榊励↓險倬鹸縲
+      !縺溘□縺励∝驛ｨ逧縺ｫ縺ｯ蛻ｻ縺ｿ蟷1e-5遘偵→縺励※縺繧九
+      r(0, :) = [zero, x] ! 髢句ｧ狗せ
+      do i = 1, max_i
+        call rktype%solver(rktype%time+1e-2,x) !遨榊螳溯｡
+        r(i, :) = [rktype%time, x]
+      end do
+    end subroutine
+end module
+