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)

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)

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)

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)

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)

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)

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)

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)

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)