Mercurial > op > rk_sakamoti

comparison rkmain.f90 @ 1:954faa9f6837

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
separate output from calcuration.
author "uncorrelated zombie" <uncorrelated@yahoo.co.jp>
date Sun, 21 Jul 2024 02:08:58 +0900
parents 5bd0f2a15d2a
children 206e90e91191
comparison
equal deleted inserted replaced
0:5bd0f2a15d2a 1:954faa9f6837
1 program rkmain 1 program rkmain
2 use rungekutta !ルンゲクッタ解法をまとめたモジュール 2 use rungekutta !ルンゲクッタ解法をまとめたモジュール
3 implicit none 3 implicit none
4 real(kind=dpkind),dimension(3) :: x !状態変数 4 real(kind=dpkind),dimension(3) :: x !状態変数
5 type(rk) :: rktype !ルンゲクッタ法ソルバに関連する型 5 type(rk) :: rktype !ルンゲクッタ法ソルバに関連する型
6 real(kind=dpkind),allocatable,dimension(:,:) :: r
7 real(kind=dpkind) :: max_t
6 8
7 !微分方程式の右辺を計算する関数、解法、時間刻みの指定
8 call rktype%setup(rhside,"rk4",1d-5)
9 !初期値の指定 9 !初期値の指定
10 x=[0.5e0,0.5e0,0.5e0] 10 x = [0.5e0,0.5e0,0.5e0]
11 !時刻0?50秒まで、0.01秒毎に標準出力へ書き出し。 11 max_t = 5e1
12 !ただし、内部的には刻み幅1e-5秒としている。 12 call calc(x, max_t, r)
13 do 13 print *, r;
14 if(rktype%time>5e1)exit 14
15 call rktype%solver(rktype%time+1e-2,x) !積分実行
16 print*,rktype%time,x
17 enddo
18 !======================================== 15 !========================================
19 contains 16 contains
20 !微分方程式の右辺を定義 17 !微分方程式の右辺を定義
21 function rhside(self,time,wk) 18 function rhside(self,time,wk)
22 class(rk),intent(inout) :: self !parameter 19 class(rk),intent(inout) :: self !parameter
26 real(kind=dpkind),parameter :: a=10e0,b=28e0,c=8e0/3e0 23 real(kind=dpkind),parameter :: a=10e0,b=28e0,c=8e0/3e0
27 rhside(1)=a*(wk(2)-wk(1)) 24 rhside(1)=a*(wk(2)-wk(1))
28 rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3) 25 rhside(2)=b*wk(1)-wk(2)-wk(1)*wk(3)
29 rhside(3)=wk(1)*wk(2)-c*wk(3) 26 rhside(3)=wk(1)*wk(2)-c*wk(3)
30 end function 27 end function
28 !Rから呼び出すサブルーチン
29 subroutine calc(x, max_t, r)
30 implicit none
31 real(kind=dpkind),dimension(3),intent(in) :: x
32 real(kind=dpkind),intent(in) :: max_t
33 real(kind=dpkind),allocatable,dimension(:,:),intent(out) :: r
34 real(kind=dpkind) :: zero = 0
35 integer i, max_i
36 !微分方程式の右辺を計算する関数、解法、時間刻みの指定
37 call rktype%setup(rhside,"rk4",1d-5)
38 max_i = int(max_t / 0.01)
39 allocate(r(0:max_i, 4))
40 !時刻max_t秒まで、0.01秒毎に配列に記録。
41 !ただし、内部的には刻み幅1e-5秒としている。
42 r(0, :) = [zero, x] ! 開始点
43 do i = 1, max_i
44 call rktype%solver(rktype%time+1e-2,x) !積分実行
45 r(i, :) = [rktype%time, x]
46 end do
47 end subroutine
31 end program 48 end program
49