An array with real subscripted variable

Question

I have problem using the results of a subroutine in the main program. I wrote this code:

Program RK4
implicit none

real k1,k2,k3,k4,h,t,R
integer i,n
real a
read*,n,h 


t=0
R=0

Do i=1,n

  call Scale_Factor(h,n,t,a)

  k1=h*(1/a(t))

  k2=h*(1/a(t+h/2.0))

  k3=h*(1/a(t+h/2.0))

  k4=h*(1/a(t+h))

  t=t+h

  R=R+(k1+2*k2+2*k3+k4)*(1/6.0)

  write(*,*)t,R

End Do

end program

!-----------------------------------------

SUBROUTINE Scale_Factor(h,n,t,a)
  implicit none
  real t,a,k1,k2,k3,k4,h,g
  integer i,n

  t=0
  a=0.001


Do i=1,n

   k1=h*g(a)

   k2=h*g(a+k1/2.0)

   k3=h*g(a+k2/2.0)

   k4=h*g(a+k3)

   t=t+h

   a=a+(k1+2*k2+2*k3+k4)*(1/6.0)

   write(*,*)t,a

END DO
END SUBROUTINE

!-------------------------
FUNCTION g(a)
  implicit none
  real a,g
  g=sqrt((1.0/a)+(1.0/a**2)) 
END FUNCTION

The subroutine solves a differential equation and produces a for each t. I need to call the result of the subroutine in the main program and use a(t) in the main program. I wanted to define a(t) as an array but since t is real, I could not do it. Can anyone help me?

Looks like you didn't specify any sizes for the variables like `a` — albert, Jun 11 '20 at 12:45
You declare `real a` so a floating point scalar, but you use `a(t)` so an array element (where `t` is also a real which is not good as index). So you have to declare your array probably as `allocatable` or `pointer`. That said actually why invent the wheel, there are on the net a lot of RK4 methods in Fortran available, just google for it. — albert, Jun 11 '20 at 13:09
You don't declare `a` to be an array, but you don't use it consistent with being a function. — francescalus, Jun 11 '20 at 13:14
@albert, `real a` in this example's program _doesn't_ declare `a` to be a "floating point scalar". It declares it to be a function (with obvious result type). — francescalus, Jun 11 '20 at 13:28
@francescalus you are probably right but I see nowhere the function `a` probably should be the `function g` as used in `Scale_Factor` but from the code in the program there is no real evidence for it (so where `a(..)` is used, probaby `g(...)` should be used). — albert, Jun 11 '20 at 13:34
No, a(t) should not be g(t). Because I want to use the results of the subroutine in the main program. Actually the function in the main program is 1/a and I need to get a from the subroutine. — Elham Q, Jun 11 '20 at 13:56
@albert is probably right that `a` wants to be declared as an array (and used correctly, with elements of `a` passed to the subroutine in the loop); it's just that it's currently without a `dimension` attribute being specified `a` is a function not an array (or scalar). — francescalus, Jun 11 '20 at 14:10
@francescalus I tried to define a(t) as an array but since t is real I couldn't — Elham Q, Jun 12 '20 at 05:59
You need to solve the coupled ODE system as a coupled ODE system, not as two strangely entwined scalar RK4 methods. — Lutz Lehmann, Jun 12 '20 at 06:51
You have to consider your arrays to be indexed by the timestep, not by the time itself. If you want to have an array with an element fr every time you have, you index the array using the integer timestep number `i`, not using the time `t`. — Vladimir F Героям слава, Jun 12 '20 at 07:43
@Vladimir F I modified the code based on what you said and it works now. But I have another problem. Actually this code is for solving two differential equations. One in the subroutine and one in the main program. I solved also these two differential equations with maple and compared its result with the result of my code. Unfortunately, they produce different results. I also checked the subroutine separately and this happened for it, too. Is there anything wrong with the algorithm in the subroutine? — Elham Q, Jun 12 '20 at 16:49
Who knows, you would have to provide much more detail. You can ask a new question, but read [ask] before doing that. Be aware that there is also a dedicated computational science stackexchange https://scicomp.stackexchange.com If you want to ask there, read their rules first. — Vladimir F Героям слава, Jun 12 '20 at 17:31

score 0 · Answer 1 · answered Jun 13 '20 at 12:54

You have a system

R'(t)=f(a(t))
a'(t)=g(a(t))

that is semi-coupled. To integrate both functions together, use a coupled RK4 method of the principal form

rk4step(R,a,h)
    k1a = h*g(a)
    k1R = h*f(a)
    k2a = h*g(a+k1a/2)
    k2R = h*f(a+k1a/2)
    k3a = h*g(a+k2a/2)
    k3R = h*f(a+k2a/2)
    k4a = h*g(a+k3a)
    k4R = h*f(a+k3a)
    a += (k1a+2*k2a+2*k3a+k4a)/6
    R += (k1R+2*k2R+2*k3R+k4R)/6
    return R,a

The better approach is to use vector-valued states and functions to avoid repetitions of similar steps.

An array with real subscripted variable

1 Answers1