В настоящее время я работаю над библиотекой временных интеграторов.Мне интересно, что лучше для вычислительной эффективности: вычисление коэффициентов для каждого временного шага или использование их в качестве параметров из внешнего репозитория.Первый случай будет таким:
Subroutine getRKF45(this, dt, y, z)
Implicit none
Class(diffEc), Intent(InOut) :: this
Real(8), Intent(In) :: dt
Real(8), dimension(:), Intent(InOut) :: y, z
Real(8), dimension(6,this%n) :: k
Call field(this%deriv, this%x, this%t)
k(1,:) = this%deriv(:)
Call field(this%deriv, this%x+k(1,:)*dt*(1.d0/5.d0), this%t+dt*(1.d0/5.d0))
k(2,:) = this%deriv(:)
Call field(this%deriv, this%x+k(1,:)*dt*(3.d0/40.d0)+k(2,:)*dt*(9.d0/40.d0) &
, this%t + dt*(3.d0/10.d0))
k(3,:) = this%deriv(:)
Call field(this%deriv, this%x+k(1,:)*dt*(3.d0/10.d0)-k(2,:)*dt*(9.d0/10.d0) &
+k(3,:)*dt*(6.d0/5.d0), this%t + dt*(3.d0/5.d0))
k(4,:) = this%deriv(:)
Call field(this%deriv, this%x-k(1,:)*dt*(11.d0/54.d0)+k(2,:)*dt*(5.d0/2.d0) &
-k(3,:)*dt*(70.d0/27.d0)+k(4,:)*dt*(35.d0/27.d0), this%t + dt)
k(5,:) = this%deriv(:)
Call field(this%deriv, this%x+k(1,:)*dt*(1631.d0/55296.d0)+k(2,:)*dt*(175.d0/512.d0) &
+k(3,:)*dt*(575.d0/13824.d0)+k(4,:)*dt*(44275.d0/110592.d0) &
+k(5,:)*dt*(253.d0/4096.d0), this%t + dt*(7.d0/8.d0))
k(6,:) = this%deriv(:)
y(:) = this%x(:) + dt*(k(1,:)*(37.d0/378.d0)+k(3,:)*(250.d0/621.d0)+k(4,:)*(125.d0/594.d0) &
+k(6,:)*(512.d0/1771.d0))
z(:) = this%x(:) + dt*(k(1,:)*(2825.d0/27648.d0)+k(3,:)*(18575.d0/48384.d0) &
+k(4,:)*(13525.d0/55296.d0)+k(5,:)*(277.d0/14336.d0)+k(6,:)*(1.d0/4.d0))
End Subroutine getRKF45
И с использованием предварительно вычисленных коэффициентов:
Subroutine getRKF45(this, dt, y, z)
Use IntegratorUtilities, only: rk45Coef
Implicit none
Class(diffEc), Intent(InOut) :: this
Real(8), Intent(In) :: dt
Real(8), dimension(:), Intent(InOut) :: y, z
Real(8), dimension(6,this%n) :: k
Call field(this%deriv, this%x, this%t)
k(1,:) = this%deriv(:)
Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(1), this%t+dt*rk45Coef(1))
k(2,:) = this%deriv(:)
Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(2)+k(2,:)*dt*rk45Coef(3) &
, this%t + dt*rk45Coef(4))
k(3,:) = this%deriv(:)
Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(4)-k(2,:)*dt*rk45Coef(5) &
+k(3,:)*dt*rk45Coef(6), this%t + dt*rk45Coef(7))
k(4,:) = this%deriv(:)
Call field(this%deriv, this%x-k(1,:)*dt*rk45Coef(8)+k(2,:)*dt*rk45Coef(9) &
-k(3,:)*dt*rk45Coef(10)+k(4,:)*dt*rk45Coef(11), this%t + dt)
k(5,:) = this%deriv(:)
Call field(this%deriv, this%x+k(1,:)*dt*rk45Coef(12)+k(2,:)*dt*rk45Coef(13) &
+k(3,:)*dt*rk45Coef(14)+k(4,:)*dt*rk45Coef(15) &
+k(5,:)*dt*rk45Coef(16), this%t + dt*rk45Coef(17))
k(6,:) = this%deriv(:)
y(:) = this%x(:) + dt*(k(1,:)*rk45Coef(18)+k(3,:)*rk45Coef(19)+k(4,:)*rk45Coef(20) &
+k(6,:)*rk45Coef(21))
z(:) = this%x(:) + dt*(k(1,:)*rk45Coef(22)+k(3,:)*rk45Coef(23) &
+k(4,:)*rk45Coef(24)+k(5,:)*rk45Coef(25)+k(6,:)*rk45Coef(26))
End Subroutine getRKF45
Я провел некоторый тест, и разница во времени между двумя схемами не имеет значения,Компилятор предварительно обрабатывает эти коэффициенты умножения?