-1

The code below works however I'm having trouble storing the final result as I'd like to. I solve for the variable x which depends on t, so I want to store x(t) as a 1D array. For each value of t, I have a different value of x.

I am having trouble doing so since I do not know how many iterations the DO WHILE loop will necessarily take, so I do not know how to write the data to a variable-size array file. As of now, I just write the results for x and t to the screen.

How can I store the solution x(t) in a single 1D array? I want to have an array such that, given an input value of time t, it gives me the value of x at this time.

program RK4

implicit none

real :: a,b,x,dt,t,f,xn,k1,k2,k3,k4,phi
integer :: n=100.

write(*,*) 'enter an initial time, final time, x(0) = '
read(*,*) a,b,x

t = a !start time here... initial time
dt = (b-a)/real(n)

do while (t < b) 

k1 = f(x,t) 
k2 = f(x+k1*dt/2., t+dt/2.) 
k3 = f(x+k2*dt/2., t+dt/2.)
k4 = f(x+k3*dt, t+dt)
phi = (1./6.)*(k1+2.*k2+2.*k3+k4)

xn = x + dt*phi
t = t + dt 
x = xn 

write(*,*) 'the results for position and time are given by', x,t !This prints out the results, but I want to store this in a single variable as x(t).

end do


end program RK4

function f(x,t)
real, intent(in) :: x,t
f = -2*x
end function f
Jeff Faraci
  • 403
  • 13
  • 28
  • Will you need the contents of the 1D array for any calculations prior to the end of the do-loop, or is this just for storage? – Matt P Aug 07 '20 at 02:27
  • I do not need the contents prior to the end of the loop. After the DO WHILE loop is over, I will need to use the data for other calculations in the code. So I’d like to have it stored nicely before I do that. Thanks for your help. – Jeff Faraci Aug 07 '20 at 02:31
  • The easiest thing might be to use a scratch file. Keep track of how many items are stored there, then once done, read it into an appropriately sized allocatable array... – Matt P Aug 07 '20 at 02:39
  • Ok that sounds reasonable. Can you provide a basic code for this? If so I’ll be able to apply it to my real code. This code is just a sample, my actual code is a bit more complicated and the function I’m trying to store is a 5D array opposed to 1D array. Thanks a lot for your help on my problems! – Jeff Faraci Aug 07 '20 at 02:43

2 Answers2

1

You want to store x(t) as a single variable -- this is not really possible1, as t is not an integer. The best you can do is to store two arrays, t_array for all the t values, and x_array with the corresponding x values.

Your values for t go from a to b in steps of dt -- and you have calculated dt as 1/100th of the distance between a and b -- so you will have exactly 101 values for x and t. (Absent of rounding errors, which you might want to protect against.) You can just make your array that large. If you want to be more flexible with the array length, you can use allocatable arrays:

real, allocatable :: t_array(:), x_array(:)
integer :: i

read (*,*) a, b, x, n

allocate(x_array(0:n), t_array(0:n))
dt = (b-a)/real(n)
t_array(:) = [(a + i*dt, i=0, n)]
x_array(0) = x

do i = 1, n
    t = t_array(i-1)    ! Makes it easier to type
    k1 = f(x, t)
    ...
    x = xn
    x_array(i) = xn
    print *, "Step: ", i, x_array(i), t_array(i)
end do

If you don't know how many values for x and t you need (for example if you want to keep going until phi is sufficiently small or something), then you can use the move_alloc command to grow your array as you populate it. But that's more convoluted, see here: How to increase array size on-the-fly in Fortran? , a specific example of this you can see in my answer to the question Fortran array input

1You can create custom types, but I think that would be overkill.

chw21
  • 7,970
  • 1
  • 16
  • 31
1

Consider that the solution to a function is computed N times, each time with a different set of independent variables. Each solution can be stored in an array.

Suppose that you do not know the number of times the function will be called - for example, perhaps the step size where your function is evaluated is dependent on convergence properties. In that case, N may exceed your initial expectations. If so, you could move your data into a new larger array (as pointed out by chw21).

Alternatively, you could store the output from each solution in a temporary "scratch" file, keeping track of how many times results are stored there. Once the solution has been obtained for every iteration and N is known, we can read the output data into an appropriately sized allocatable array. Then, the output from each iteration can be accessed (by an index number 1<=index<=N) as needed. Note that the scratch file will be automatically deleted once it is closed.

For example:

integer :: tmpfile                 !! scratch file logical unit number
integer :: n                       !! number of iterations
integer :: i                       !! index counter
real :: x,y,z                      !! output results 
real, allocatable :: results(:,:)  !! 2D array to store multiple values from each iteration

n = 0
open(newunit=tmpfile, status="scratch")
do while (...condition...)
  ...(get function output from the current iteration)...
  n = n + 1                        !! update the number of output iterations
  write(tmpfile, *) x,y,z          !! write your output values to the scratch file
enddo

allocate(results(n,3))             !! n is known. allocate storage array
rewind(tmpfile)                    !! important! return to the beginning of the scratch file
do i = 1,n
  read(tmpfile,*) results(i,:)     !! simplest, but list-directed 'read' is not only option
enddo
close(tmpfile)

  
Matt P
  • 2,287
  • 1
  • 11
  • 26