Skip to content

Commit

Permalink
finish explicit single pass runge_kutta routines
Browse files Browse the repository at this point in the history
  • Loading branch information
@cbcoutinho
cbcoutinho committed Nov 6, 2016
1 parent 7dcd42f commit a1ea03f
Show file tree
Hide file tree
Showing 5 changed files with 80 additions and 89 deletions.
1 change: 1 addition & 0 deletions .gitignore
View file
Original file line number Diff line number Diff line change
@@ -1,4 +1,5 @@
.ipynb_checkpoints
data.out

main

Expand Down
100 changes: 50 additions & 50 deletions data.out
View file
Original file line number Diff line number Diff line change
@@ -1,50 +1,50 @@
0.0000000000000000 0.0000000000000000 0.10000000000000001
0.25645654315018718 2.5153696175667262E-002 9.4537193820639354E-002
0.51291308630037435 4.7559197013290616E-002 8.3045725841833534E-002
0.76936962945056153 6.5850199825126660E-002 6.6546202878540808E-002
1.0258261726007487 7.8991760728469598E-002 4.6347153600479876E-002
1.2822827157509360 8.6334616144943518E-002 2.3945950932185794E-002
1.5387392589011231 8.7641615143106630E-002 9.2148300542800946E-004
1.7951958020513104 8.3085710610873290E-002 -2.1173961421678902E-002
2.0516523452014974 7.3220865353094688E-002 -4.0916396161173180E-002
2.3081088883516849 5.8928965419590605E-002 -5.7098966751420781E-002
2.5645654315018720 4.1347289639118261E-002 -6.8802573788480925E-002
2.8210219746520586 2.1782176973920590E-002 -7.5444394148701577E-002
3.0774785178022461 1.6152051784527020E-003 -7.6802035739631236E-002
3.3339350609524336 -1.7791581076519389E-002 -7.3012773188573521E-002
3.5903916041026207 -3.5185072623348163E-002 -6.4548986654270363E-002
3.8468481472528078 -4.9499436289354541E-002 -5.2172454355741332E-002
4.1033046904029948 -5.9918678681022845E-002 -3.6871436482379553E-002
4.3597612335531828 -6.5919966507826569E-002 -1.9785458983902292E-002
4.6162177767033699 -6.7295660743776589E-002 -2.1233096134345014E-003
4.8726743198535560 -6.4153520079256038E-002 1.4920028718232500E-002
5.1291308630037440 -5.6895998926982556E-002 3.0241957994182264E-002
5.3855874061539311 -4.6180910550360860E-002 4.2901345150320520E-002
5.6420439493041172 -3.2866862900744753E-002 5.2173933747414407E-002
5.8985004924543052 -1.7947737105522329E-002 5.7591003707171552E-002
6.1549570356044923 -2.4810209168472513E-003 5.8959438059830949E-002
6.4114135787546802 1.2484990363560112E-002 5.6362666697660699E-002
6.6678701219048673 2.5980253476120530E-002 5.0143246914837195E-002
6.9243266650550543 3.7173882565386250E-002 4.0869024498789805E-002
7.1807832082052414 4.5423215663880662E-002 2.9285823425828895E-002
7.4372397513554285 5.0308300478805260E-002 1.6260377992735425E-002
7.6936962945056155 5.1650141608825494E-002 2.7177080053804419E-003
7.9501528376558026 4.9512198495985232E-002 -1.0422674810060324E-002
8.2066093808059897 4.4185755101908863E-002 -2.2307452266055100E-002
8.4630659239561776 3.6160824174269121E-002 -3.2203189976465665E-002
8.7195224671063656 2.6085135871216303E-002 -3.9539775971892940E-002
8.9759790102565518 1.4714440340439441E-002 -4.3941170470850756E-002
9.2324355534067397 2.8577904679937625E-003 -4.5241975112269676E-002
9.4888920965569259 -8.6783540498962912E-003 -4.3489333631828130E-002
9.7453486397071121 -1.9143447239164089E-002 -3.8930668816066223E-002
10.001805182857300 -2.7890179973588922E-002 -3.1988677278227233E-002
10.258261726007488 -3.4412928191736886E-002 -2.3225786705773147E-002
10.514718269157674 -3.8375160446358435E-002 -1.3300883591905004E-002
10.771174812307862 -3.9624463657542165E-002 -2.9215108355350131E-003
11.027631355458048 -3.8194723967629109E-002 7.2051028385751309E-003
11.284087898608234 -3.4295870549376155E-002 1.6418886857432348E-002
11.540544441758424 -2.8292396698194212E-002 2.4148633970392522E-002
11.797000984908610 -2.0672563883589522E-002 2.9946024407905020E-002
12.053457528058798 -1.2010729790068928E-002 3.3510045047075963E-002
12.309914071208985 -2.9255919814400108E-003 3.4700598715412055E-002
12.566370614359172 5.9627106101172167E-003 3.3540866782862486E-002
0.0000000000000000 0.0000000000000000 1.0000000000000000
0.25645654315018718 0.25645654315018718 0.96711502073772815
0.51291308630037435 0.49604595009403868 0.86954150481199244
0.76936962945056153 0.70273309776175763 0.71373242085379796
1.0258261726007487 0.86266508380132967 0.51004084398600680
1.2822827157509360 0.96509967212429215 0.27203205612055975
1.5387392589011231 1.0031267901591758 1.5580161887986943E-002
1.7951958020513104 0.97413462092687542 -0.24219062035816774
2.0516523452014974 0.87998885483857892 -0.48404938427599709
2.3081088883516849 0.72691280779078304 -0.69381033003475923
2.5645654315018720 0.52507609633849839 -0.85741593757723678
2.8210219746520586 0.28791905240635735 -0.96387903280863385
3.0774785178022461 3.1257755569713219E-002 -1.0060206156905953
3.3339350609524336 -0.22777072451181496 -0.98095390454623421
3.5903916041026207 -0.47185253630929719 -0.89028196309887297
3.8468481472528078 -0.68465431012321498 -0.73999538886626126
4.1033046904029948 -0.85191612670870787 -0.54007657822214150
4.3597612335531828 -0.96240705483592903 -0.30383670623765557
4.6162177767033699 -1.0086792301596335 -4.7029456067338271E-002
4.8726743198535560 -0.98756984632281264 0.21319949513428366
5.1291308630037440 -0.90041722688293835 0.45945718306537603
5.3855874061539311 -0.75297622415503151 0.67526583252764549
5.6420439493041172 -0.55503827552123630 0.84616540914953298
5.8985004924543052 -0.31978119757720025 0.96068247467341616
6.1549570356044923 -6.2891893006697952E-002 1.0111004319112287
6.4114135787546802 0.19847962713620726 0.99397945264842291
6.6678701219048673 0.44686516310237934 0.91039105989986535
6.9243266650550543 0.66564575561748029 0.76585137379101975
7.1807832082052414 0.84016360463724693 0.56995715779766631
7.4372397513554285 0.95870408435433074 0.33574867475723769
7.6936962945056155 1.0132820649171810 7.8841651173033545E-002
7.9501528376558026 1.0001797625416364 -0.18361387049553873
8.2066093808059897 0.92019989329019869 -0.43407837660223059
8.4630659239561776 0.77861689896255259 -0.65579500182984107
8.7195224671063656 0.58482917920246291 -0.83391059514062005
8.9759790102565518 0.35173525524633975 -0.95647073214445688
9.2324355534067397 9.4875271081657886E-002 -1.0152220196174415
9.4888920965569259 -0.16860502992140547 -1.0061678486077590
9.7453486397071121 -0.42109878529173206 -0.92984017664055063
10.001805182857300 -0.64571455785344056 -0.79126886291383836
10.258261726007488 -0.82740632529434543 -0.59965027939719040
10.514718269157674 -0.95398132319881568 -0.36773702642901113
10.771174812307862 -1.0169182337551410 -0.11098924956351380
11.027631355458048 -1.0119408179965692 0.15345596450195492
11.284087898608234 -0.93930837900015640 0.40792841234760907
11.540544441758424 -0.80380333195244602 0.63540547479745246
11.797000984908610 -0.61441638448498670 0.82065080282073088
12.053457528058798 -0.38375004639796728 0.95123482020826033
12.309914071208985 -0.12718004036557395 1.0183686932051073
12.566370614359172 0.13816958733617818 1.0174958133575667
17 changes: 8 additions & 9 deletions src/main.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,26 +6,25 @@ program main
use lib_constants, only: pi => pi_dp
implicit none

integer :: ii, ios
integer, parameter :: N = 50
real(wp) :: dt
real(wp), dimension(N) :: t
real(wp), dimension(2) :: y0
real(wp), dimension(N,2) :: y
integer :: ii, ios
integer, parameter :: N = 50, num_eq = 2
real(wp) :: dt
real(wp), dimension(N) :: t
real(wp), dimension(num_eq) :: y0
real(wp), dimension(N,num_eq) :: y

call linspace(0._wp, 4._wp*pi, t)
dt = t(2) - t(1)
y(1,:) = [0._wp, 0.1_wp]
y(1,:) = [0._wp, 1._wp]

do ii = 1, N-1
y0 = y(ii, :)
call rk(t(ii), dt, y0, mysub)
call rk(num_eq, t(ii), dt, y0, mysub)
y(ii+1, :) = y0
end do

open(unit=21, file='data.out', iostat=ios, status="replace", action="write")
if ( ios /= 0 ) stop "Error opening file 21"

do ii = 1, N
write(21,*) t(ii), y(ii, :)
end do
Expand Down
7 changes: 4 additions & 3 deletions src/misc.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,10 +4,11 @@ module misc

contains

subroutine mysub(t, y, dy)
subroutine mysub(N, t, y, dy)
integer, intent(in) :: N
real(wp), intent(in) :: t
real(wp), intent(in), dimension(2) :: y
real(wp), intent(out), dimension(2) :: dy
real(wp), intent(in), dimension(N) :: y
real(wp), intent(out), dimension(N) :: dy

! dy(1) = dcos(t)
! dy(2) = dsin(t)
Expand Down
44 changes: 17 additions & 27 deletions src/runge_kutta.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -4,66 +4,56 @@ module runge_kutta

contains

subroutine rk(t, dt, y, dysub)
subroutine rk(n, t, dt, y, dysub)
! Dummy arguments
integer, intent(in) :: n
real(wp), intent(in) :: t, dt
real(wp), intent(inout), dimension(2) :: y
real(wp), dimension(2) :: dy
real(wp), intent(inout), dimension(n) :: y
real(wp), dimension(n) :: dy

interface
subroutine dysub(t, y, dy)
subroutine dysub(n, t, y, dy)
import wp
integer, intent(in) :: n
real(wp), intent(in) :: t
real(wp), intent(in), dimension(2) :: y
real(wp), intent(out), dimension(2) :: dy
real(wp), intent(in), dimension(n) :: y
real(wp), intent(out), dimension(n) :: dy
end subroutine
end interface

! Local arguments
integer :: ii, jj, n, m
integer :: ii, jj, m
real(wp) :: dummy_t
real(wp), dimension(size(y)) :: dummy_y
real(wp), dimension(n) :: dummy_y
real(wp), dimension(:), allocatable :: b, c
real(wp), dimension(:,:), allocatable :: a, k


call classic_rk4(a, b, c)
call heun(a, b, c)

n = size(b)
m = size(y)
! print *, c
! print *, b
!
! do ii = 1, n
! print *, a(ii, :)
! end do
m = size(b)

allocate(k(n, m))
allocate(k(m, n))
k = 0._wp

call dysub(t, y, dy)
call dysub(n, t, y, dy)
k(1, :) = dy

! print *, t
! print *, y
! print *, dy
! print *,

dummy_y = y

do ii = 2, N
do ii = 2, n

dummy_t = t + dt*c(ii)
do jj = 1, ii-1
dummy_y = dummy_y + dt*a(ii,jj)*k(jj,:)
end do

call dysub(dummy_t, dummy_y, dy)
call dysub(n, dummy_t, dummy_y, dy)
k(ii, :) = dy

end do

y = y + [( dt*dot_product(b,k(:,jj)), jj = 1, m )]
y = y + [( dt*dot_product(b,k(:,jj)), jj = 1, n )]

deallocate(a, b, c, k)

Expand Down

0 comments on commit a1ea03f

Please sign in to comment.