使用 Fortran 求解微分方程使用 RK4 比使用 C++ 快 5 倍

Question

为什么我的 Fortran 代码执行速度比使用 RK4 求解这个二阶微分方程（行星和太阳之间的万有引力）的 C++ 代码快 5 倍？请问如何优化我的 C++ 代码？

我已经尝试将 pow() 更改为 x*x，但没有任何改进。删除写操作将 Fortran 上的执行时间除以 2，但只使 C++ 代码快了大约 15%。

代码如下：

C++（编译：c++ -Wall -Wextra equadiff.cpp -o equadiffcpp.x）：

#include<iostream>
#include<fstream>
#include<cmath>
#include <chrono> //pour mesurer le temps d'execution

using namespace std;
using namespace std::chrono;

void deriv(double t, double X[], double Xderiv[], int n){
double radius;
radius=sqrt(pow(X[0],2)+pow(X[1],2));
Xderiv[0]=X[2];
Xderiv[1]=X[3];
Xderiv[2]=-X[0]/pow(radius,3);
Xderiv[3]=-X[1]/pow(radius,3);
}



void rk4(double t, double X[], double dt, int n, void deriv(double, double[], double[], int)){

int i;
double ddt;
double Xp[n], k1[n], k2[n], k3[n], k4[n];
ddt=0.5*dt;

deriv(t,X,k1,n);

for(i=0;i<n;i++){
Xp[i]=X[i]+ddt*k1[i];
deriv(t+ddt,Xp,k2,n);
}

for(i=0;i<n;i++){
Xp[i]=X[i]+ddt*k2[i];
deriv(t+ddt,Xp,k3,n);
}

for(i=0;i<n;i++){
Xp[i]=X[i]+dt*k3[i];
deriv(t+dt,Xp,k4,n);
}

for(i=0;i<n;i++){
X[i] = X[i] + dt*(k1[i]+2*k2[i]+2*k3[i]+k4[i])/6;
}
}


int main(void){
double dt, t, tmax;
double X[4];
double Xderiv[4];

dt=0.01;
tmax=1000.0;

X[0]=1.0;
X[1]=0.0;
X[2]=0.0;
X[3]=-0.5;

ofstream fichier ("equadiffrk4cpp.out");

auto start = high_resolution_clock::now();//on commence a compter le temps de mesure

for(t=0.0;t<=tmax;t+=dt){
    rk4(t,X,dt,4,deriv);
    if((int)(round(t/dt))%10==0){//on n'ecrit qu'une valeur sur 10
    fichier <<t<<" "<<X[0]<<" "<<X[1]<<endl;
    }
}

auto stop = high_resolution_clock::now();//on arrete de compter le temps d'execution

fichier.close();

auto duration = duration_cast<microseconds>(stop - start); 
  
cout << "Time taken by function: "
         << duration.count() << " microseconds" << endl; 
  
return 0; 

}

Fortran90（编译：gfortran equadiff.f90 -o equadifff90.x）：

program  meca_planet
implicit none
real (8) :: dt ,t
integer  :: i
real(8),  dimension (4) :: X, Xderiv
external  :: euler, deriv_planet, rk4

real :: start, finish!pour compter temps execution du programme

t=0.
dt=0.01

!Initialization
X=(/1.,0.,0.,-0.5/)

open(11,file='equadiffrk4f90.out')

call cpu_time(start)!on commence a compter

do i=1,100000!tmax=0.01*100000=1000
    t=t+dt
    call  rk4(t,X,dt,4,deriv_planet)
    if (mod(nint(t/dt),10).eq.0) then
        write(11,*) t, X(1),X(2)
    endif
enddo

call cpu_time(finish)!on arrete de compter

close (11)

print '("Time = ",f6.3," seconds.")',finish-start

end program meca_planet

subroutine deriv_planet(t,X,Xderiv,n)
implicit none
integer, intent(in) :: n
real(8), intent (in) :: t!pourquoi on definit t dans deriv_planet mais ensuite on ne l'utilise pas?
real(8) :: radius
real(8), dimension(n), intent(in) :: X
real(8), dimension(n), intent(out) :: Xderiv
if (n.ne.4) write (*,*) 'WARNING: dimension de n incorrecte, devrait etre 4'
radius=sqrt(X(1)**2+X(2)**2)
Xderiv(1)=X(3)
Xderiv(2)=X(4)
Xderiv(3)=-X(1)/radius**3
Xderiv(4)=-X(2)/radius**3
end subroutine deriv_planet

subroutine rk4(t,X,dt,n,deriv)
!Runge-Kutta du 4eme ordre
implicit none
integer, intent(in) :: n
real(8), intent(in) :: t, dt
real(8), dimension(n), intent (inout) :: X
real(8) :: ddt
real(8), dimension(n) :: Xp, k1, k2, k3, k4
ddt = 0.5*dt
call deriv(t,X,k1,n); Xp=X+ddt*k1
call deriv(t+ddt,Xp,k2,n); Xp=X+ddt*k2
call deriv(t+ddt,Xp,k3,n); Xp=X+dt*k3
call deriv(t+dt,Xp,k4,n); X = X + dt*(k1+2.0*k2 + 2.0*k3+k4)/6.0
end subroutine rk4

最终目标是为太阳系编写一个 N 体问题，然后可能是一个星系。我正在考虑使用 C++，但基于这个初步评估，我现在倾向于使用 Fortran90。

Answer 1

Bob__ 的评论应该抓住了罪魁祸首：在您的循环中，您在每个坐标更新后调用 deriv。但只有最后一次通话才算数。正如使用完全设置的矢量完成的那样，您总体上会得到正确的结果。改为

for(i=0;i<n;i++) {
  Xp[i]=X[i]+ddt*k1[i];
}
deriv(t+ddt,Xp,k2,n);

等并且时间应该减少大约 3 倍，将 3*4+1=13 deriv 调用替换为 4.

您可以将根计算减少一个平方根。只需计算

r3 = pow(X[0]*X[0]+X[1]*X[1], 1.5);

或许也可以使用

用乘法代替除法

ir3 = pow(X[0]*X[0]+X[1]*X[1], -1.5);

用RK4模拟一个大太阳系是不切实际的想法。你会得到太多不同的尺度，频率大不相同的过程。此外，明显偏心轨道上的物体在远点移动非常缓慢，在最近点移动非常快，这很难用固定的步长捕获。最好使用具有自适应步长的求解器，甚至可能对通过“密集输出”插值连接的内部系统、外部系统和彗星使用单独的求解器。查看 "Moving stars around" 教程，最后他们给出了一些如何组织大规模天体模拟的初步想法。

使用 Fortran 求解微分方程使用 RK4 比使用 C++ 快 5 倍

Solving differential equation using RK4 is 5 times faster with Fortran than with C++

c++

performance

differential-equations