• Complex RK4 integration: problems with external functions

    From Khilav Jayesh Majmudar@21:1/5 to All on Sun Nov 7 08:43:23 2021
    Hello All,

    I am trying to perform a 4th order Runge-Kutta integration of two coupled complex-valued differential equations in Fortran 90 using the gcc compiler. Following are what I think the relevant snippets of my code:

    !!!!! RK4 integration using functions fb and fe which are defined later. The increment is over an index k. The increment is from k to k+2 because of the way the system is constructed and is not important for this question.

    do k=1,994,2
    k1b(k) = dz1*fb(f1(k))
    k1e(k) = dz1*fe(f2(k))

    k2b(k) = dz1*fb(f1(k)+k1b(k)/2.0)
    k2e(k) = dz1*fe(f2(k)+k1e(k)/2.0)

    k3b(k) = dz1*fb(f1(k)+k2b(k)/2.0)
    k3e(k) = dz1*fe(f2(k)+k2e(k)/2.0)

    k4b(k) = dz1*fb(f1(k)+k3b(k))
    k4e(k) = dz1*fe(f2(k)+k3e(k))

    f2(k+2) = f2(k)+(1.0/6.0)*(k1b(k)+2*(k2b(k)+k3b(k))+k4b(k))
    f1(k+2) = f1(k)+(1.0/6.0)*(k1e(k)+2*(k2e(k)+k3e(k))+k4e(k))
    end do

    !!!!! The external functions 1are given below. array1 and array2 are globally defined arrays. constant is a globally defined real number.

    complex(8) function fb(f1)

    complex(8) :: ij = (0,1)
    complex(8),intent(in) :: f1

    fb = (ij*constant-array1)*f1/(array2)**2

    end function fb

    complex(8) function fe(f2)

    complex(8) :: ij = (0,1)
    complex(8),intent(in) :: f2

    fe = (ij*constant)*f2

    end function fe

    The given code compiles but doesn't increment the functions of interest (f1 and f2). The first value itself (k=1) is nonsensical and thereafter I get NaNs. The initial values given to them are sensible.

    I suspect this is because f1 and f2 are arrays and so they need to be defined as such inside the external functions fe and fb. However, when I do so, gcc gives me an incompatible rank error between fe and f2, since I defined both fe and fb as complex
    scalars.

    But when I define fe and fb as complex arrays, gcc shoots out an error which says that "array index must be of type integer" inside the various k functions, which makes sense.

    So I need help in marrying these two issues: to get fe and fb to read f1 and f2 as arrays without being arrays themselves.

    Any help would be greatly appreciated. Thank you!

    Khilav

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Thomas Koenig@21:1/5 to Khilav Jayesh Majmudar on Sun Nov 7 18:14:07 2021
    Khilav Jayesh Majmudar <majmu008@umn.edu> schrieb:

    Hello All,

    I am trying to perform a 4th order Runge-Kutta integration of
    two coupled complex-valued differential equations in Fortran 90
    using the gcc compiler. Following are what I think the relevant
    snippets of my code:

    You left out some all-important details, i.e. the declaration of
    your variables and functions. Without these, it is only possible
    to give some general advice. If you have not declared them,
    the type inferred from Fortran's implicit type rules is almost
    certainly wrong.

    The best way is to put your functions fb and fe into a module,
    like this:

    module x
    implicit none
    integer, parameter :: wp = selected_real_kind(15) ! double precision
    contains
    complex (kind=wp) function fb(f1)
    fb = ...
    end function fb
    complex (kind=wp) function fe(f2)
    fe = ...
    end function fe
    end module x

    !!!!! RK4 integration using functions fb and fe which are defined
    later. The increment is over an index k. The increment is from
    k to k+2 because of the way the system is constructed and is not
    important for this question.



    do k=1,994,2
    k1b(k) = dz1*fb(f1(k))
    k1e(k) = dz1*fe(f2(k))

    k2b(k) = dz1*fb(f1(k)+k1b(k)/2.0)
    k2e(k) = dz1*fe(f2(k)+k1e(k)/2.0)

    k3b(k) = dz1*fb(f1(k)+k2b(k)/2.0)
    k3e(k) = dz1*fe(f2(k)+k2e(k)/2.0)

    k4b(k) = dz1*fb(f1(k)+k3b(k))
    k4e(k) = dz1*fe(f2(k)+k3e(k))

    f2(k+2) = f2(k)+(1.0/6.0)*(k1b(k)+2*(k2b(k)+k3b(k))+k4b(k))
    f1(k+2) = f1(k)+(1.0/6.0)*(k1e(k)+2*(k2e(k)+k3e(k))+k4e(k))
    end do

    !!!!! The external functions 1are given below. array1 and array2 are globally defined arrays. constant is a globally defined real number.

    How did you globally define an array? Put it in a module?

    complex(8) function fb(f1)

    complex(8) :: ij = (0,1)
    complex(8),intent(in) :: f1

    fb = (ij*constant-array1)*f1/(array2)**2

    You have to know what you want to write.

    Is fb an array-valued function or not? Assigning

    fb = (ij*constant-array1)*f1/(array2)**2

    if array1 and array2 gives you an array, which you cannot assign
    to a scalar function. Do you want to return a scalar or
    an array?

    A minor point, 1.0/6.0 will give you default REAL precision,
    which is less than double precision that you are otherwise
    using.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From pehache@21:1/5 to All on Sun Nov 7 19:17:15 2021
    Le 07/11/2021 à 17:43, Khilav Jayesh Majmudar a écrit :
    Hello All,

    I am trying to perform a 4th order Runge-Kutta integration of two coupled complex-valued differential equations in Fortran 90 using the gcc compiler. Following are what I think the relevant snippets of my code:

    !!!!! RK4 integration using functions fb and fe which are defined later. The increment is over an index k. The increment is from k to k+2 because of the way the system is constructed and is not important for this question.

    do k=1,994,2
    k1b(k) = dz1*fb(f1(k))
    k1e(k) = dz1*fe(f2(k))

    k2b(k) = dz1*fb(f1(k)+k1b(k)/2.0)
    k2e(k) = dz1*fe(f2(k)+k1e(k)/2.0)

    k3b(k) = dz1*fb(f1(k)+k2b(k)/2.0)
    k3e(k) = dz1*fe(f2(k)+k2e(k)/2.0)

    k4b(k) = dz1*fb(f1(k)+k3b(k))
    k4e(k) = dz1*fe(f2(k)+k3e(k))

    f2(k+2) = f2(k)+(1.0/6.0)*(k1b(k)+2*(k2b(k)+k3b(k))+k4b(k))
    f1(k+2) = f1(k)+(1.0/6.0)*(k1e(k)+2*(k2e(k)+k3e(k))+k4e(k))
    end do

    !!!!! The external functions 1are given below. array1 and array2 are globally defined arrays. constant is a globally defined real number.

    complex(8) function fb(f1)

    complex(8) :: ij = (0,1)
    complex(8),intent(in) :: f1

    fb = (ij*constant-array1)*f1/(array2)**2

    This is inconsistent, as the RHS expression is an array, and the LHS
    (fb) a scalar

    end function fb

    complex(8) function fe(f2)

    complex(8) :: ij = (0,1)
    complex(8),intent(in) :: f2

    fe = (ij*constant)*f2

    end function fe

    The given code compiles but doesn't increment the functions of interest (f1 and f2). The first value itself (k=1) is nonsensical and thereafter I get NaNs. The initial values given to them are sensible.

    I suspect this is because f1 and f2 are arrays and so they need to be defined as such inside the external functions fe and fb.

    No, because you don't pass the full f1 and f2 arrays to fb() and fe(),
    but only scalar elements (f1(k) and f2(k))

    However, when I do so, gcc gives me an incompatible rank error between fe and f2, since I defined both fe and fb as complex scalars.

    indeed...


    But when I define fe and fb as complex arrays,

    Defining fe() and fb() as arrays would be inconsistent with the main
    code, where you assign the results of fe() and fb() to scalars.



    So I need help in marrying these two issues: to get fe and fb to read f1 and f2 as arrays without being arrays themselves.


    Looks like you should first write mathematically what you want to do...


    --
    "...sois ouvert aux idées des autres pour peu qu'elles aillent dans le
    même sens que les tiennes.", ST sur fr.bio.medecine

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Khilav Jayesh Majmudar@21:1/5 to All on Sun Nov 7 17:55:25 2021
    You left out some all-important details, i.e. the declaration of
    your variables and functions. Without these, it is only possible
    to give some general advice. If you have not declared them,
    the type inferred from Fortran's implicit type rules is almost
    certainly wrong.
    Here are the declarations:

    real(8) :: constant
    complex(8),dimension(996) :: cEx2,cBy2,k1b,k1e,k2b,k2e,k3b,k3e,k4b,k4e complex(8) :: fb,fe
    real(8),parameter :: dz1=4.0E4

    The best way is to put your functions fb and fe into a module,
    like this:

    module x
    implicit none
    integer, parameter :: wp = selected_real_kind(15) ! double precision contains
    complex (kind=wp) function fb(f1)
    fb = ...
    end function fb
    complex (kind=wp) function fe(f2)
    fe = ...
    end function fe
    end module x
    Thank you. I will try out this approach.

    How did you globally define an array? Put it in a module?
    array1 and array2 have been computed earlier in the main code. Currently there are no modules being used.

    You have to know what you want to write.

    Is fb an array-valued function or not? Assigning

    fb = (ij*constant-array1)*f1/(array2)**2

    if array1 and array2 gives you an array, which you cannot assign
    to a scalar function. Do you want to return a scalar or
    an array?
    The way I have written the code right now, I want to return a scalar. I would prefer it to be an array as I would be able to keep store all values of k1b,k1e,etc. The problem is, when I tried to define fb and fe as arrays as I mentioned in my original
    post, I wasn't able to continue as the arguments of fb and fe become complex numbers during the k computations inside the loop.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From pehache@21:1/5 to All on Mon Nov 8 11:34:40 2021
    Le 08/11/2021 à 02:55, Khilav Jayesh Majmudar a écrit :


    if array1 and array2 gives you an array, which you cannot assign
    to a scalar function. Do you want to return a scalar or
    an array?
    The way I have written the code right now, I want to return a scalar. I would prefer it to be an array as I would be able to keep store all values of k1b,k1e,etc. The problem is, when I tried to define fb and fe as arrays as I mentioned in my original post, I wasn't able to continue as the arguments of fb
    and fe become complex numbers during the k computations inside the loop.

    As f1 and f2 are updated during the iterations, you probably wouldn't want
    to return full arrays in fb() and fe()

    Would you really want to do in your computation is quite unclear, but I
    feel that you need to modify fb() to pass the proper elements of array1
    and array2 :

    ====================================
    complex(8) function fb(f1,v1,v2)

    complex(8) :: ij = (0,1)
    complex(8),intent(in) :: f1
    complex(8),intent(in) :: v1, v2

    fb = (ij*constant-v1)*f1/(v2)**2

    end function fb
    ====================================

    And the calls to fb() would be something like (3 arguments) :

    k1b(k) = dz1*fb(f1(k),array1(k),array2(k))

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Thomas Koenig@21:1/5 to Khilav Jayesh Majmudar on Mon Nov 8 20:14:07 2021
    Khilav Jayesh Majmudar <majmu008@umn.edu> schrieb:
    You left out some all-important details, i.e. the declaration of
    your variables and functions. Without these, it is only possible
    to give some general advice. If you have not declared them,
    the type inferred from Fortran's implicit type rules is almost
    certainly wrong.
    Here are the declarations:

    real(8) :: constant
    complex(8),dimension(996) :: cEx2,cBy2,k1b,k1e,k2b,k2e,k3b,k3e,k4b,k4e complex(8) :: fb,fe
    real(8),parameter :: dz1=4.0E4

    The best way is to put your functions fb and fe into a module,
    like this:

    module x
    implicit none
    integer, parameter :: wp = selected_real_kind(15) ! double precision
    contains
    complex (kind=wp) function fb(f1)
    fb = ...
    end function fb
    complex (kind=wp) function fe(f2)
    fe = ...
    end function fe
    end module x
    Thank you. I will try out this approach.

    How did you globally define an array? Put it in a module?
    array1 and array2 have been computed earlier in the main
    code. Currently there are no modules being used.

    Generally, variables declared in different parts have nothing to
    do with each other, even if they share the same name.

    So,

    program main
    implicit none
    real, dimension(10) :: asdf
    asdf(1) = 42.
    call foo
    print *,asdf(1)
    end program main

    subroutine foo
    implicit none
    real, dimension(10) :: asdf
    asdf(1) = 43.
    end subroutine foo

    will happily print 42 because the arrays are totally different
    entities.

    If you want to share as a global variable, either put them
    in COMMON (but you don't want that), or use a module:

    module x
    implicit none
    real, dimension(10) :: asdf
    contains
    subroutine foo
    asdf(1) = 43.
    end subroutine foo
    end module x

    program main
    use x
    implicit none
    asdf(1) = 42.
    call foo
    print *,asdf(1)
    end program main

    This will print 43, as expected.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Khilav Jayesh Majmudar@21:1/5 to All on Tue Nov 9 12:26:45 2021
    Would you really want to do in your computation is quite unclear, but I
    feel that you need to modify fb() to pass the proper elements of array1
    and array2 :

    ====================================
    complex(8) function fb(f1,v1,v2)
    complex(8) :: ij = (0,1)
    complex(8),intent(in) :: f1
    complex(8),intent(in) :: v1, v2

    fb = (ij*constant-v1)*f1/(v2)**2

    end function fb
    ====================================

    And the calls to fb() would be something like (3 arguments) :

    k1b(k) = dz1*fb(f1(k),array1(k),array2(k))

    Thank you. I realised the silly mistake of not declaring the input of v1 and v2 inside my function. The code compiles now. However, the values of the k1b,k1e,k2b, etc. are not getting incremented over the k loop.

    I am trying to integrate a set of coupled differential equations using the 4th order Runge-Kutta method (https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method). In that link, the k1,k2,k3,k4 you see are my k1b,k2b, etc.
    The f functions given there are defined as fb and fe in my case. It is a fairly standard way of solving ordinary differential equations so I think I am making a stupid error somewhere.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From Robin Vowels@21:1/5 to majm...@umn.edu on Tue Nov 9 21:42:10 2021
    On Wednesday, November 10, 2021 at 7:26:47 AM UTC+11, majm...@umn.edu wrote:
    Would you really want to do in your computation is quite unclear, but I feel that you need to modify fb() to pass the proper elements of array1 and array2 :

    ====================================
    complex(8) function fb(f1,v1,v2)
    complex(8) :: ij = (0,1)
    complex(8),intent(in) :: f1
    complex(8),intent(in) :: v1, v2

    fb = (ij*constant-v1)*f1/(v2)**2

    end function fb
    ====================================

    And the calls to fb() would be something like (3 arguments) :

    k1b(k) = dz1*fb(f1(k),array1(k),array2(k))
    Thank you. I realised the silly mistake of not declaring the input of v1 and v2 inside my function. The code compiles now. However, the values of the k1b,k1e,k2b, etc. are not getting incremented over the k loop.

    I am trying to integrate a set of coupled differential equations using the 4th order Runge-Kutta method (https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method). In that link, the k1,k2,k3,k4 you see are my k1b,k2b,
    etc. The f functions given there are defined as fb and fe in my case. It is a fairly standard way of solving ordinary differential equations so I think I am making a stupid error somewhere.
    .
    This really isn't going anywhere, with only a few lines of code posted.
    .
    Better to post the entire code.

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)
  • From gah4@21:1/5 to majm...@umn.edu on Wed Nov 10 11:14:01 2021
    On Sunday, November 7, 2021 at 8:43:25 AM UTC-8, majm...@umn.edu wrote:

    I am trying to perform a 4th order Runge-Kutta integration of two coupled complex-valued differential equations in Fortran 90 using the gcc compiler. Following are what I think the relevant snippets of my code:

    !!!!! RK4 integration using functions fb and fe which are defined later. The increment is over an index k. The increment is from k to k+2 because of the way the system is constructed and is not important for this question.

    (snip)

    !!!!! The external functions 1are given below. array1 and array2 are globally defined arrays. constant is a globally defined real number.

    There aren't global variables in Fortran. There is COMMON, which is sometimes called global, but needs to be declared.
    For internal procedures, you get access to host variables. I presume this is what you mean, but the details aren't shown.

    complex(8) function fb(f1)

    complex(8) :: ij = (0,1)
    complex(8),intent(in) :: f1
    fb = (ij*constant-array1)*f1/(array2)**2
    end function fb

    Assuming array1 and array2 are arrays, this doesn't make sense.

    (snip)

    The given code compiles but doesn't increment the functions of interest (f1 and f2).
    The first value itself (k=1) is nonsensical and thereafter I get NaNs. The initial values given to them are sensible.

    You have to initialize the k=1 elements with the initial value. (That is, this is an initial value problem.)

    If they are being changed, then something is very wrong.

    I suspect this is because f1 and f2 are arrays and so they need to be defined as such inside the external functions fe and fb. However, when I do so, gcc gives me an incompatible rank error between fe and f2, since I defined both fe and fb as complex
    scalars.

    When I write these, I use scalar variables for the k's. They are not needed between iterations,
    though for debugging purposes, it is sometimes nice. Then I usually just print them out
    inside the loop, but could also be assigned to array elements, but only for debugging.
    (Partly because it might be faster as scalars.)

    But when I define fe and fb as complex arrays, gcc shoots out an error which says that "array index must be of type integer" inside the various k functions, which makes sense.

    You can write functions that return arrays. That isn't usual, but could be used for
    solving mutliple equations at the same time. In this case, you are solving two equations
    (the one with b and the one with e), which could be done with two element arrays.
    As they are independent, I would probably give them separate loops (it sometimes
    helps the optimization), but otherwise it works both ways.

    So I need help in marrying these two issues: to get fe and fb to read
    f1 and f2 as arrays without being arrays themselves.

    It isn't so obvious what should depend on what.

    Can you write the differential equations in math form, so it is more
    obvious what is variable and what isn't?

    --- SoupGate-Win32 v1.05
    * Origin: fsxNet Usenet Gateway (21:1/5)