Search code examples
fortrannumerical-methodscomplex-numbers

How to compute a complex function in a cleverer way with Fortran?

I want to solve this equation by using Fortran:

eq1

Where ψ (psi) is a complex Fortran variable. Now, I am solving this by defining two new complex variables:
ir=(1.0,0.0) and ii=(0.0,1.0).
I use these to select only the real or imaginary part of the equation. In this way I solve my equation separately for the real and imaginary part. The code is here:

do i = 1,nn    
  mod2 = (abs(psi(i)))**2
  psi(i) = ir*(-beta*imag(der2(i)) + alpha*mod2*imag(psi(i))) + ii*(beta*real(der2(i)) - alpha*mod2*real(psi(i)))
end do

Where psi and der2 are complex arrays with nn elements.
I want to solve this equation in a better way without splitting it in two equations. I tried to solve it in this way:

mod2 = abs(psi)**2
psi = -ii*(-beta*der2+alpha*mod2*psi)

but it doesn't work because I obtain completely different values with respect to the first method I used. For me it makes sense that it doesn't work because in the second method I am not evaluating the real part. Is this right?
As an example, the 10° element of my psi array becomes:\

  1. (-6.39094355774850412E-003,-6.04041029332164168E-003) (with 1° method)\
  2. (-1.75266632213431602E-004,-6.21567692553507290E-003) (2° method)\

Any suggestion? Thank you!

Solution

  • The problem is mod2 = abs(psi)**2 should have been mod2 = abs(dat)**2

    But I think that the correct calculation for mod2 is sum( real( dat*conjg(dat) ) )

    I get the same result with both methods with some arbitrary data:

    eq1=
 (-595.000000000000,-120.200000000000)
 (-713.800000000000,-1.40000000000000)
 (-832.600000000000,117.400000000000)
 (-951.400000000000,236.200000000000)
 (-1070.20000000000,355.000000000000)
 (-1189.00000000000,473.800000000000)
 (-1307.80000000000,592.600000000000)
 (-1426.60000000000,711.400000000000)
 (-1545.40000000000,830.200000000000)
 (-1664.20000000000,949.000000000000)
 eq2=
 (-595.000000000000,-120.200000000000)
 (-713.800000000000,-1.40000000000000)
 (-832.600000000000,117.400000000000)
 (-951.400000000000,236.200000000000)
 (-1070.20000000000,355.000000000000)
 (-1189.00000000000,473.800000000000)
 (-1307.80000000000,592.600000000000)
 (-1426.60000000000,711.400000000000)
 (-1545.40000000000,830.200000000000)
 (-1664.20000000000,949.000000000000)

    Test code

    program Console1
use,intrinsic :: iso_fortran_env
implicit none
! Variables
integer, parameter :: sp=real32, wp=real64
complex(wp), parameter :: ir = (1d0,0d0), ii = (0d0,1d0)
real(wp), parameter :: alpha = 0.1d0, beta = 0.2d0
integer, parameter :: n=10
complex(wp) :: dat(n), psi(n), der2(n)
integer :: i

! Body of Console1
    
    dat = [ ( (2d0-i)*ir - (4d0+i)*ii, i=1,n ) ]
    der2 = [ ( -(5d0+i)*ir + (1d0-i)*ii, i=1,n ) ]

    psi = eq1(dat, der2)
    
    print *, "eq1="
    do i=1,n
    print *, psi(i)
    end do

    psi = eq2(dat, der2)
    
    print *, "eq2="
    do i=1,n
    print *, psi(i)
    end do
    

contains

function eq1(dat, der2) result(psi)
complex(wp), intent(in) :: dat(:), der2(size(psi))
complex(wp) :: psi(size(dat))
real(wp) :: mod2        
integer :: i
    mod2 = sum( real( dat*conjg(dat) ) )
    do i=1, size(psi)
        psi(i) = ir*(-beta*imag(der2(i)) + alpha*mod2*imag(dat(i))) + ii*(beta*real(der2(i)) - alpha*mod2*real(dat(i)))
    end do
end function

function eq2(dat, der2) result(psi)
complex(wp), intent(in) :: dat(:), der2(size(dat))
complex(wp) :: psi(size(dat))
real(wp) :: mod2        
    mod2 = sum( real( dat*conjg(dat) ) )
    psi = -ii*(-beta*der2+alpha*mod2*dat)    
end function

end program Console1

    Also with your definition of |ψ| we have also consistent results.

    eq1=
 (-13.0000000000000,-3.80000000000000)
 (-21.4000000000000,-1.40000000000000)
 (-34.6000000000000,3.40000000000000)
 (-53.8000000000000,11.8000000000000)
 (-80.2000000000000,25.0000000000000)
 (-115.000000000000,44.2000000000000)
 (-159.400000000000,70.6000000000000)
 (-214.600000000000,105.400000000000)
 (-281.800000000000,149.800000000000)
 (-362.200000000000,205.000000000000)
 eq2=
 (-13.0000000000000,-3.80000000000000)
 (-21.4000000000000,-1.40000000000000)
 (-34.6000000000000,3.40000000000000)
 (-53.8000000000000,11.8000000000000)
 (-80.2000000000000,25.0000000000000)
 (-115.000000000000,44.2000000000000)
 (-159.400000000000,70.6000000000000)
 (-214.600000000000,105.400000000000)
 (-281.800000000000,149.800000000000)
 (-362.200000000000,205.000000000000)

    with code

    function eq1(dat, der2) result(psi)
complex(wp), intent(in) :: dat(:), der2(size(psi))
complex(wp) :: psi(size(dat))
real(wp) :: mod2(size(dat))
integer :: i
    mod2 = abs(dat)**2 
    do i=1, size(psi)
        psi(i) = ir*(-beta*imag(der2(i)) + alpha*mod2(i)*imag(dat(i))) + ii*(beta*real(der2(i)) - alpha*mod2(i)*real(dat(i)))
    end do
end function

function eq2(dat, der2) result(psi)
complex(wp), intent(in) :: dat(:), der2(size(dat))
complex(wp) :: psi(size(dat))
real(wp) :: mod2(size(dat))     
    mod2 = abs(dat)**2 
    psi = -ii*(-beta*der2+alpha*mod2*dat)    
end function

    Going back to the original problem, I can replicate the issue

     eq1=
 (-13.0000000000000,-3.80000000000000)
 (-21.4000000000000,-1.40000000000000)
 (-34.6000000000000,3.40000000000000)
 (-53.8000000000000,11.8000000000000)
 (-80.2000000000000,25.0000000000000)
 (-115.000000000000,44.2000000000000)
 (-159.400000000000,70.6000000000000)
 (-214.600000000000,105.400000000000)
 (-281.800000000000,149.800000000000)
 (-362.200000000000,205.000000000000)
 eq2=
 (0.000000000000000E+000,-1.20000000000000)
 (0.200000000000000,-1.40000000000000)
 (0.400000000000000,-1.60000000000000)
 (0.600000000000000,-1.80000000000000)
 (0.800000000000000,-2.00000000000000)
 (-1952.64000000000,779.256000000000)
 (NaN,NaN)
 (1.40000000000000,-2.60000000000000)
 (NaN,NaN)
 (1.80000000000000,-3.00000000000000)

    with mod2 = abs(psi)**2 instead of mod2 = abs(dat)**2