如何使用 Fortran 以更聪明的方式计算复杂函数？

Question

我想用 Fortran 解这个方程：

其中 ψ (psi) 是一个复杂的 Fortran 变量。现在，我通过定义两个新的复杂变量来解决这个问题：
ir=(1.0,0.0) 和 ii=(0.0,1.0).
我只用这些来 select 等式的实部或虚部。通过这种方式，我分别为实部和虚部求解方程。代码在这里：

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

其中 psi 和 der2 是具有 nn 个元素的复数数组。
我想以更好的方式解决这个方程，而不是将它分成两个方程。我试过这样解决：

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

但它不起作用，因为我获得了与我使用的第一种方法完全不同的值。对我来说，它不起作用是有道理的，因为在第二种方法中我没有评估实部。这样对吗？
例如，我的 psi 数组的 10° 元素变为：\

(-6.39094355774850412E-003,-6.04041029332164168E-003)（1°法）\
(-1.75266632213431602E-004,-6.21567692553507290E-003) (2°法)\

有什么建议吗？谢谢！

Answer 1

问题是 mod2 = abs(psi)**2 应该是 mod2 = abs(dat)**2

但我认为 mod2 的正确计算是 sum( real( dat*conjg(dat) ) )

我用这两种方法和一些任意数据得到了相同的结果：

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)

测试代码

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

此外，根据您对 |ψ| 的定义，我们也得到了一致的结果。

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)

带代码

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

回到最初的问题，我可以复制这个问题

 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)

用mod2 = abs(psi)**2代替mod2 = abs(dat)**2

如何使用 Fortran 以更聪明的方式计算复杂函数？

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

fortran

numerical-methods

complex-numbers