Варианты выбора
Подходит к концу моя студенческая лицензия на Mathematica 9 и я задумался о приобретении собственной персональной лицензии на один из программных продуктов. Выбор профессиональных инструментов не особенно велик. Собственно, всего три варианта:
Maple Personal edition - $299.00
Mathematica Home Edition - $295.00
Matlab + Image Processing Toolbox + Symbolic Math Toolbox +
Partial Differential Equation Toolbox - $135 + $39 + $39 +$39 = $252
В этом отношении Matlab выглядит привлекательнее по цене, и остаётся запас для приобретения ещё одного тулбокса, правда, несколько смущает тот факт, что за ту же цену MapleSoft и Wolfram предлагают гораздо более полный набор инструментов, хотя что мне до них, если большую их часть я всё равно не буду использовать.
Поэтому следующим вопросом стал вопрос удобства использования и скорости работы.
Для начала я написал коротенький тест на скорость умножения матриц. Здесь matlab откровенно порадовал, оставив всех далеко позади.
Matlab 2015a:
% File: TestSpeed.m
N = 5000;
tic % Запускаем "таймер"
A = rand (N, N); % Создаём матрицу nxn и заполняем случайными числами
AM=A'*A; % Произведение исходной и транспонированной матриц
toc % Останавливаем "таймер"
Время вычислений: 4.3 с.
Mathematica 9:
(* File: TestSpeed.nb *)
n = 5000;
Start = AbsoluteTime[]; (* Запускаем "таймер" *)
A = Table[RandomReal[], {i, n}, {j, n}];(* Создаём матрицу nxn
и заполняем случайными числами *)
AT = Transpose[A]; (* Транспонируем матрицу *)
AM = AT.A; (* Вычисляем произведение исходной и транспонированной матриц *)
AbsoluteTime[] - Start (* Останавливаем "таймер" *)
Время вычислений: 21.9 с.
Gfortran 4.9.1:
! File: TestSpeed.f90
program main
integer, parameter :: n = 5000
real a(n,n),at(n,n),am (n,n),ss ! Три матрицы nxn и переменная "таймера"
ss = secnds(0.) ! Запускаем "таймер"
call random_number(a) ! Заполняем матрицу случайными числами
at = transpose(a) ! Транспонируем матрицу
am = matmul(at,a) ! Произведение исходной и транспонированной матриц
ss = secnds(0.) - ss ! Останавливаем "таймер"
write(6,"('Processing time', f8.3, ' sec')"), ss
end program main
Время вычислений: 82.3 с.
Maple 18:
restart : with(LinearAlgebra) :
start := time() :
A := RandomMatrix(1000, 1000) :
MatrixMatrixMultiply(Transpose(A), A) :
time()-start
Время вычислений: 91.2 с.
Увы, дождаться же перемножения матриц 5000х5000 мне не удалось.
Результат Maple не лезет ни в какие ворота, с ним ещё предстоит разбираться отдельно, а вот фортран откровенно огорчил, но тут, видимо, проблема не в самом фортране, а в неоптимальной реализации библиотечных процедур работы с матрицами. Очевидно, что переписав их самому, можно достичь такого же, если не лучшего результата, что и в Matlab.
Однако мне в реальной жизни не нужно перемножать такие матрицы, а если и нужно, то не больше, чем 4х4. Основные численные расчеты, с которыми я работаю, относятся к решению дифференциальных уравнений. Поэтому для следующего теста я выбрал классическую задачу двух тел, в первую очередь, потому что для неё известно точное аналитическое решение.
И вот, что получилось...
На этот раз начну с фортрана, поскольку на нем у меня давно написаны все рабочие алгоритмы. Задача решалась модифицированным методом Рунге-Кутты-Фельберга (Runge-Kutta-Fehlberg), более известного как RK45, с переменным шагом интегрирования.
Я привожу программу целиком:
Gfortran 4.9.1:
! File: orbit.f90
Program RK45F
real*8 RK45 ! Объявление имени функции
real*8 y(8), T, Tmax, eps, M1, M2, V, dt
real*8 prec
common M1, M2, dt
M1 = 50. ! Масса первого тела
M2 = 10. ! Масса второго тела
V = 4. ! Начальная скорость первого тела
y(1) = 0. ! Начальная координата x первого тела
y(2) = 0. ! Начальная координата y первого тела
y(3) = 5. ! Начальная координата x второго тела
y(4) = 0. ! Начальная координата y второго тела
y(5) = 0. ! Начальный компонент скорости Vx первого тела
y(6) = -V*M2/M1 ! Начальный компонент скорости Vy первого тела
y(7) = 0. ! Начальный компонент скорости Vx второго тела
y(8) = V ! Начальный компонент скорости Vy второго тела
T=0. ! Начальный момент времени
Tmax = 401*1000. ! 1000 периодов обращения двойной системы
prec=0.05 ! Точность интегрирования
dt = 1. ! Начальный шаг интегрирования
open (1, file = 'orbit.dat') ! Открываем файл для записи данных
TOTAL = SECNDS(0.) ! Запускаем "таймер"
do while (T<=Tmax) ! Главный цикл
eps = RK45(y) ! Возвращаем максимальное из ускорений
dt = prec/eps ! Чем больше ускорение, тем меньше шаг
write (1,*), (y(i), i=1,4),eps ! Пишем данные в файл
T = T + dt ! Приращение времени
end do !
TOTAL = SECNDS(0.0) - TOTAL ! Останавливаем "таймер"
WRITE(6,1) Total
1 FORMAT('FORTRAN: Total time = ', F6.2, ' sec')
call system ('gnuplot orbit.plt')
! Вызов внешней программы gnuplot, которая в качестве параметра
! получает файл orbit.plt с инструкциями:
!
! set terminal windows
! set grid
! set linetype 1 lc rgb "red" lw 1 pt 0
! set linetype 2 lc rgb "blue" lw 1 pt 7
! set title "Orbits" textcolor lt 1
! plot 'orbit.dat' using 1:2 with lines lw 1 title "M1", 'orbit.dat' using 3:4 with lines lw 1 title "M2"
! pause 1000
!
! При запуске под Windows необходимо создать в каталоге командный файл
! gnuplot.cmd, в котором прописать полный путь к программе wgnuplot.exe.
! Например, у меня этот файл состоит из единственной строки:
!
! "C:\Program Files (x86)\Maxima-5.31.2\gnuplot\bin\wgnuplot.exe" %1 %2 %3
!
! При запуске под Linux командный файл создавать не надо, командный
! процессор сам найдет программу gnuplot (если она установлена), но
! из файла orbit.plt необходимо удалить первую и последнюю строки
! (set terminal windows и pause 1000).
close(1)
End program RK45F
Module Derivatives
! Система дифференциальных уравнений.
! Функция возвращает максимальное значение модуля вторых производных.
Implicit NONE
Contains
real*8 Function gr(f,y,FLAG)
! y(8) - входной массив
! f(8) - выходной массив
real*8 M1, M2, r12, f(8), y(8), dt
common M1, M2, dt
integer FLAG ! Флаг: вычислять или нет максимальный модуль ускорения
r12 = sqrt((y(1)-y(3))**2+(y(2)-y(4))**2) ! Расстояние между компонентами
r12=1.0/(r12**3) ! Оптимизация: замена четырёх делений
! одним делением и четырьмя умножениями.
f(1)=y(5) ! Компонент скорости Vx первого тела
f(2)=y(6) ! Компонент скорости Vy первого тела
f(3)=y(7) ! Компонент скорости Vx второго тела
f(4)=y(8) ! Компонент скорости Vy второго тела
f(5)=M2*(y(3)-y(1))*r12 ! Компонент ускорения ax первого тела
f(6)=M2*(y(4)-y(2))*r12 ! Компонент ускорения ay первого тела
f(7)=M1*(y(1)-y(3))*r12 ! Компонент ускорения ax второго тела
f(8)=M1*(y(2)-y(4))*r12 ! Компонент ускорения ay второго тела
if (FLAG.eq.1) then
gr = dmax1(abs(f(5)),abs(f(6)),abs(f(7)),abs(f(8)))
! Максимальный модуль компонента ускорения
else
gr=0
endif
End Function gr
End Module Derivatives
real*8 Function RK45 (y) ! Реализация модифицированного алгоритма RK45.
! Вычисление погрешности на каждом шаге не
! используется, вместо этого в основной программе
! производится автоматическое изменение шага.
use Derivatives
real*8 M1, M2, dt
common M1, M2, dt
real*8 k1(8), k2(8), k3(8), k4(8), k5(8), y(8), y1(8), f(8)
real*8 eps
real*8, parameter :: c21=1./4.
real*8, parameter :: c31=3./32., c32=9./32.
real*8, parameter :: c41=1932./2197., c42=-7200./2197., c43=7296./2197.
real*8, parameter :: c51=439./216.,c52=-8.,c53=3680./513.,c54=-845./4104.
real*8, parameter :: a1=25./216.,a3=1408./2565.,a4=2197./4104.,a5=-1./5.
integer i
eps = gr(f,y,0)
do i=1,8
k1(i) = dt*f(i)
y1(i) = y(i)+c21*k1(i)
end do
eps = gr(f,y1,0)
do i=1,8
k2(i) = dt*f(i)
y1(i) = y(i)+c31*k1(i)+c32*k2(i)
end do
eps = gr(f,y1,0)
do i=1,8
k3(i) = dt*f(i)
y1(i) = y(i)+c41*k1(i)+c42*k2(i)+c43*k3(i)
end do
eps = gr(f,y1,0)
do i=1,8
k4(i) = dt*f(i)
y1(i) = y(i)+c51*k1(i)+c52*k2(i)+c53*k3(i)+c54*k4(i)
end do
eps = gr(f,y1,1) ! Модуль ускорения вычисляется только
! при последнем обращении к функции
do i=1,8
k5(i) = dt*f(i)
y(i) = y(i)+a1*k1(i)+a3*k3(i)+a4*k4(i)+a5*k5(i)
end do
RK45 = eps
End function RK45
Время работы фортрановской программы: 7.6 с.
Фортран как честная числомолотилка ожидаемо занял первое место.
Результат работы приведен на следующем рисунке:
Теперь переходим к Mathematica 9
(* File orbit.nb *)
M1 = 50; M2 = 10; V = 4.0; T = 1000; Start = AbsoluteTime[];
System = {
x1'[t]==vx1[t],
vx1'[t]==M2(x2[t]-x1[t])/((x2[t]-x1[t])^2 + (y2[t]-y1[t])^2)^(3/2),
y1'[t]==vy1[t],
vy1'[t]==M2(y2[t]-y1[t])/((x2[t]-x1[t])^2 + (y2[t]-y1[t])^2)^(3/2),
x2'[t]==vx2[t],
vx2'[t]==M1(x1[t]-x2[t])/((x1[t]-x2[t])^2 + (y1[t]-y2[t])^2)^(3/2),
y2'[t]==vy2[t],
vy2'[t]==M1(y1[t]-y2[t])/((x1[t]-x2[t])^2 + (y1[t]-y2[t])^2)^(3/2)
};
Solution = NDSolve[{System,
x1[0] == 0, y1[0] == 0,
vx1[0] == 0, vy1[0] == -V*M2/M1,
x2[0] == 5, y2[0] == 0,
vx2[0] == 0, vy2[0] == V},
{x1, vx1, y1, vy1, x2, vx2, y2, vy2},
{t, 0, 401*T}, MaxSteps -> 1000000];
ParametricPlot[{{{x1[t], y1[t]} /. Solution}, {{x2[t], y2[t]} /.Solution}},
{t, 0, 401*T}, PlotPoints -> 10000]
AbsoluteTime[] - Start
Время вычислений: 16.9 с.
Результат приведен на следующем рисунке:
Для Maple 18 я задачу двух тел писать не стал, потому что до этого запускал на нём задачу трёх тел, и она считалась за точно такое же время, как и в Mathematica. Поэтому я предполагаю, что время расчёта в Maple должно составить около 17 секунд.
И наконец, Matlab 2015a:
% File: orbit.m
function orbit()
global M1;
global M2;
M1=50.0; M2=10.0; V = 4.0;
tic
[t,h]=ode45(@body2,[0,401*1000],[0,0,5,0,0,-V*M2/M1,0,V]);
x1=h(:,1); y1=h(:,2);
x2=h(:,3); y2=h(:,4);
plot(x1,y1,x2,y2);
toc
function f=body2(t,x)
global M1;
global M2;
r12=sqrt((x(1)-x(3))^2+(x(2)-x(4))^2);
r12=1.0/r12^3;
f=[x(5);x(6);x(7);x(8);...
M2*(x(3)-x(1))*r12;...
M2*(x(4)-x(2))*r12;...
M1*(x(1)-x(3))*r12;...
M1*(x(2)-x(4))*r12];
Время работы: 8.4 с.
Производительность матлаба оказалась, на первый взгляд, на высоте, но результат несколько обескураживает:
Очевидно, что функция ode45 в Matlab требует дополнительной настройки. Для того, чтобы получить такой же результат, который дают Mathematica и фортрановская программа, пришлось несколько модифицировать матлабовский сценарий (добаленный текст выделен полужирным):
% File: orbit.m
function orbit()
global M1;
global M2;
M1=50.0; M2=10.0; V = 4.0;
tic
options=odeset('RelTol',0.0000001,'AbsTol',0.0000001,...
'InitialStep',1.0, 'MaxStep', 1000000);
[t,h]=ode45(@body2,[0,401*1000],[0,0,5,0,0,-V*M2/M1,0,V],options);
x1=h(:,1); y1=h(:,2);
x2=h(:,3); y2=h(:,4);
plot(x1,y1,x2,y2);
toc
function f=body2(t,x)
global M1;
global M2;
r12=sqrt((x(1)-x(3))^2+(x(2)-x(4))^2);
r12=1.0/r12^3;
f=[x(5);x(6);x(7);x(8);...
M2*(x(3)-x(1))*r12;...
M2*(x(4)-x(2))*r12;...
M1*(x(1)-x(3))*r12;...
M1*(x(2)-x(4))*r12];
Время расчета увеличилось до 22.6 с
Это примерно соответствует времени Mathematica и Maple. Но, по крайней мере, результат стал похож на правду:
Несколько неожиданным оказался тот факт, что функция численного решения дифференциальных уравнений требует в Matlab дополнительных настроек, в то время как в Mathematica и Maple всё корректно считается с первого раза с настройками по умолчанию. Очевидно, что при численном моделировании желательно иметь достаточно точное представление о том, что именно делает та или иная функция, и с этой точки зрения мне, пожалуй, больше подойдёт фортрановская программа, полную реализацию которой я вижу перед глазами. Да и время работы ее в три раза меньше, чем время работы матлабовского скрипта.
Ради очистки совести я запустил оба приведенных здесь матлабовских скрипта на Octave 4.0.0 (в скрипт для решения задачи двух тел в начало функции orbit() необходимо добавить строку: pkg load odepkg).
Результаты получились следующие:
Умножение матриц 5000x5000: 7.0 с.
Задача двух тел (1000 периодов): 327.3 с
Последний результат, к сожалению, делает Octave совершенно неприемлемым для моих задач. Если кто-нибудь знает, как ускорить odepkg в octave, буду благодарен за совет.
С другой стороны, при запуске первого варианта скрипта для задачи двух тел Octave сразу же выдал мне предупреждение:
warning: Option "RelTol" not set, new value 0.000001 is used
warning: called from
ode45 at line 113 column 5
orbit at line 10 column 5
warning: Option "AbsTol" not set, new value 0.000001 is used
warning: Option "InitialStep" not set, new value 40100.000000 is used
warning: Option "MaxStep" not set, new value 40100.000000 is used
С подставленными им значениями по умолчанию скрипт посчитался за 160 секунд и выдал, хотя и не вполне корректный, но, по крайней мере, похожий на правду результат:
Теперь я сведу все результаты в одну табличку:
Умножение матриц Задача двух тел
GFortran 82,3 с 7,6 с
Maple N/A ~17,0 c
Mathematica 21,9 с 16,9 с
Matlab 4,3 с 22,6 с
Octave 7,0 с 327,3 с
Я, конечно, ещё поиграюсь с матлабом до истечения триального периода, но пока склоняюсь к тому, чтобы вообще отказаться от проприетарных систем компьютерной алгебры; для решения численных задач небесной механики продолжить использовать фортран, для работы с матрицами - Octave, а для символьной математики - Maxima.
Maple Personal edition - $299.00
Mathematica Home Edition - $295.00
Matlab + Image Processing Toolbox + Symbolic Math Toolbox +
Partial Differential Equation Toolbox - $135 + $39 + $39 +$39 = $252
В этом отношении Matlab выглядит привлекательнее по цене, и остаётся запас для приобретения ещё одного тулбокса, правда, несколько смущает тот факт, что за ту же цену MapleSoft и Wolfram предлагают гораздо более полный набор инструментов, хотя что мне до них, если большую их часть я всё равно не буду использовать.
Поэтому следующим вопросом стал вопрос удобства использования и скорости работы.
Для начала я написал коротенький тест на скорость умножения матриц. Здесь matlab откровенно порадовал, оставив всех далеко позади.
Matlab 2015a:
% File: TestSpeed.m
N = 5000;
tic % Запускаем "таймер"
A = rand (N, N); % Создаём матрицу nxn и заполняем случайными числами
AM=A'*A; % Произведение исходной и транспонированной матриц
toc % Останавливаем "таймер"
Время вычислений: 4.3 с.
Mathematica 9:
(* File: TestSpeed.nb *)
n = 5000;
Start = AbsoluteTime[]; (* Запускаем "таймер" *)
A = Table[RandomReal[], {i, n}, {j, n}];(* Создаём матрицу nxn
и заполняем случайными числами *)
AT = Transpose[A]; (* Транспонируем матрицу *)
AM = AT.A; (* Вычисляем произведение исходной и транспонированной матриц *)
AbsoluteTime[] - Start (* Останавливаем "таймер" *)
Время вычислений: 21.9 с.
Gfortran 4.9.1:
! File: TestSpeed.f90
program main
integer, parameter :: n = 5000
real a(n,n),at(n,n),am (n,n),ss ! Три матрицы nxn и переменная "таймера"
ss = secnds(0.) ! Запускаем "таймер"
call random_number(a) ! Заполняем матрицу случайными числами
at = transpose(a) ! Транспонируем матрицу
am = matmul(at,a) ! Произведение исходной и транспонированной матриц
ss = secnds(0.) - ss ! Останавливаем "таймер"
write(6,"('Processing time', f8.3, ' sec')"), ss
end program main
Время вычислений: 82.3 с.
Maple 18:
restart : with(LinearAlgebra) :
start := time() :
A := RandomMatrix(1000, 1000) :
MatrixMatrixMultiply(Transpose(A), A) :
time()-start
Время вычислений: 91.2 с.
Увы, дождаться же перемножения матриц 5000х5000 мне не удалось.
Результат Maple не лезет ни в какие ворота, с ним ещё предстоит разбираться отдельно, а вот фортран откровенно огорчил, но тут, видимо, проблема не в самом фортране, а в неоптимальной реализации библиотечных процедур работы с матрицами. Очевидно, что переписав их самому, можно достичь такого же, если не лучшего результата, что и в Matlab.
Однако мне в реальной жизни не нужно перемножать такие матрицы, а если и нужно, то не больше, чем 4х4. Основные численные расчеты, с которыми я работаю, относятся к решению дифференциальных уравнений. Поэтому для следующего теста я выбрал классическую задачу двух тел, в первую очередь, потому что для неё известно точное аналитическое решение.
И вот, что получилось...
На этот раз начну с фортрана, поскольку на нем у меня давно написаны все рабочие алгоритмы. Задача решалась модифицированным методом Рунге-Кутты-Фельберга (Runge-Kutta-Fehlberg), более известного как RK45, с переменным шагом интегрирования.
Я привожу программу целиком:
Gfortran 4.9.1:
! File: orbit.f90
Program RK45F
real*8 RK45 ! Объявление имени функции
real*8 y(8), T, Tmax, eps, M1, M2, V, dt
real*8 prec
common M1, M2, dt
M1 = 50. ! Масса первого тела
M2 = 10. ! Масса второго тела
V = 4. ! Начальная скорость первого тела
y(1) = 0. ! Начальная координата x первого тела
y(2) = 0. ! Начальная координата y первого тела
y(3) = 5. ! Начальная координата x второго тела
y(4) = 0. ! Начальная координата y второго тела
y(5) = 0. ! Начальный компонент скорости Vx первого тела
y(6) = -V*M2/M1 ! Начальный компонент скорости Vy первого тела
y(7) = 0. ! Начальный компонент скорости Vx второго тела
y(8) = V ! Начальный компонент скорости Vy второго тела
T=0. ! Начальный момент времени
Tmax = 401*1000. ! 1000 периодов обращения двойной системы
prec=0.05 ! Точность интегрирования
dt = 1. ! Начальный шаг интегрирования
open (1, file = 'orbit.dat') ! Открываем файл для записи данных
TOTAL = SECNDS(0.) ! Запускаем "таймер"
do while (T<=Tmax) ! Главный цикл
eps = RK45(y) ! Возвращаем максимальное из ускорений
dt = prec/eps ! Чем больше ускорение, тем меньше шаг
write (1,*), (y(i), i=1,4),eps ! Пишем данные в файл
T = T + dt ! Приращение времени
end do !
TOTAL = SECNDS(0.0) - TOTAL ! Останавливаем "таймер"
WRITE(6,1) Total
1 FORMAT('FORTRAN: Total time = ', F6.2, ' sec')
call system ('gnuplot orbit.plt')
! Вызов внешней программы gnuplot, которая в качестве параметра
! получает файл orbit.plt с инструкциями:
!
! set terminal windows
! set grid
! set linetype 1 lc rgb "red" lw 1 pt 0
! set linetype 2 lc rgb "blue" lw 1 pt 7
! set title "Orbits" textcolor lt 1
! plot 'orbit.dat' using 1:2 with lines lw 1 title "M1", 'orbit.dat' using 3:4 with lines lw 1 title "M2"
! pause 1000
!
! При запуске под Windows необходимо создать в каталоге командный файл
! gnuplot.cmd, в котором прописать полный путь к программе wgnuplot.exe.
! Например, у меня этот файл состоит из единственной строки:
!
! "C:\Program Files (x86)\Maxima-5.31.2\gnuplot\bin\wgnuplot.exe" %1 %2 %3
!
! При запуске под Linux командный файл создавать не надо, командный
! процессор сам найдет программу gnuplot (если она установлена), но
! из файла orbit.plt необходимо удалить первую и последнюю строки
! (set terminal windows и pause 1000).
close(1)
End program RK45F
Module Derivatives
! Система дифференциальных уравнений.
! Функция возвращает максимальное значение модуля вторых производных.
Implicit NONE
Contains
real*8 Function gr(f,y,FLAG)
! y(8) - входной массив
! f(8) - выходной массив
real*8 M1, M2, r12, f(8), y(8), dt
common M1, M2, dt
integer FLAG ! Флаг: вычислять или нет максимальный модуль ускорения
r12 = sqrt((y(1)-y(3))**2+(y(2)-y(4))**2) ! Расстояние между компонентами
r12=1.0/(r12**3) ! Оптимизация: замена четырёх делений
! одним делением и четырьмя умножениями.
f(1)=y(5) ! Компонент скорости Vx первого тела
f(2)=y(6) ! Компонент скорости Vy первого тела
f(3)=y(7) ! Компонент скорости Vx второго тела
f(4)=y(8) ! Компонент скорости Vy второго тела
f(5)=M2*(y(3)-y(1))*r12 ! Компонент ускорения ax первого тела
f(6)=M2*(y(4)-y(2))*r12 ! Компонент ускорения ay первого тела
f(7)=M1*(y(1)-y(3))*r12 ! Компонент ускорения ax второго тела
f(8)=M1*(y(2)-y(4))*r12 ! Компонент ускорения ay второго тела
if (FLAG.eq.1) then
gr = dmax1(abs(f(5)),abs(f(6)),abs(f(7)),abs(f(8)))
! Максимальный модуль компонента ускорения
else
gr=0
endif
End Function gr
End Module Derivatives
real*8 Function RK45 (y) ! Реализация модифицированного алгоритма RK45.
! Вычисление погрешности на каждом шаге не
! используется, вместо этого в основной программе
! производится автоматическое изменение шага.
use Derivatives
real*8 M1, M2, dt
common M1, M2, dt
real*8 k1(8), k2(8), k3(8), k4(8), k5(8), y(8), y1(8), f(8)
real*8 eps
real*8, parameter :: c21=1./4.
real*8, parameter :: c31=3./32., c32=9./32.
real*8, parameter :: c41=1932./2197., c42=-7200./2197., c43=7296./2197.
real*8, parameter :: c51=439./216.,c52=-8.,c53=3680./513.,c54=-845./4104.
real*8, parameter :: a1=25./216.,a3=1408./2565.,a4=2197./4104.,a5=-1./5.
integer i
eps = gr(f,y,0)
do i=1,8
k1(i) = dt*f(i)
y1(i) = y(i)+c21*k1(i)
end do
eps = gr(f,y1,0)
do i=1,8
k2(i) = dt*f(i)
y1(i) = y(i)+c31*k1(i)+c32*k2(i)
end do
eps = gr(f,y1,0)
do i=1,8
k3(i) = dt*f(i)
y1(i) = y(i)+c41*k1(i)+c42*k2(i)+c43*k3(i)
end do
eps = gr(f,y1,0)
do i=1,8
k4(i) = dt*f(i)
y1(i) = y(i)+c51*k1(i)+c52*k2(i)+c53*k3(i)+c54*k4(i)
end do
eps = gr(f,y1,1) ! Модуль ускорения вычисляется только
! при последнем обращении к функции
do i=1,8
k5(i) = dt*f(i)
y(i) = y(i)+a1*k1(i)+a3*k3(i)+a4*k4(i)+a5*k5(i)
end do
RK45 = eps
End function RK45
Время работы фортрановской программы: 7.6 с.
Фортран как честная числомолотилка ожидаемо занял первое место.
Результат работы приведен на следующем рисунке:
Теперь переходим к Mathematica 9
(* File orbit.nb *)
M1 = 50; M2 = 10; V = 4.0; T = 1000; Start = AbsoluteTime[];
System = {
x1'[t]==vx1[t],
vx1'[t]==M2(x2[t]-x1[t])/((x2[t]-x1[t])^2 + (y2[t]-y1[t])^2)^(3/2),
y1'[t]==vy1[t],
vy1'[t]==M2(y2[t]-y1[t])/((x2[t]-x1[t])^2 + (y2[t]-y1[t])^2)^(3/2),
x2'[t]==vx2[t],
vx2'[t]==M1(x1[t]-x2[t])/((x1[t]-x2[t])^2 + (y1[t]-y2[t])^2)^(3/2),
y2'[t]==vy2[t],
vy2'[t]==M1(y1[t]-y2[t])/((x1[t]-x2[t])^2 + (y1[t]-y2[t])^2)^(3/2)
};
Solution = NDSolve[{System,
x1[0] == 0, y1[0] == 0,
vx1[0] == 0, vy1[0] == -V*M2/M1,
x2[0] == 5, y2[0] == 0,
vx2[0] == 0, vy2[0] == V},
{x1, vx1, y1, vy1, x2, vx2, y2, vy2},
{t, 0, 401*T}, MaxSteps -> 1000000];
ParametricPlot[{{{x1[t], y1[t]} /. Solution}, {{x2[t], y2[t]} /.Solution}},
{t, 0, 401*T}, PlotPoints -> 10000]
AbsoluteTime[] - Start
Время вычислений: 16.9 с.
Результат приведен на следующем рисунке:
Для Maple 18 я задачу двух тел писать не стал, потому что до этого запускал на нём задачу трёх тел, и она считалась за точно такое же время, как и в Mathematica. Поэтому я предполагаю, что время расчёта в Maple должно составить около 17 секунд.
И наконец, Matlab 2015a:
% File: orbit.m
function orbit()
global M1;
global M2;
M1=50.0; M2=10.0; V = 4.0;
tic
[t,h]=ode45(@body2,[0,401*1000],[0,0,5,0,0,-V*M2/M1,0,V]);
x1=h(:,1); y1=h(:,2);
x2=h(:,3); y2=h(:,4);
plot(x1,y1,x2,y2);
toc
function f=body2(t,x)
global M1;
global M2;
r12=sqrt((x(1)-x(3))^2+(x(2)-x(4))^2);
r12=1.0/r12^3;
f=[x(5);x(6);x(7);x(8);...
M2*(x(3)-x(1))*r12;...
M2*(x(4)-x(2))*r12;...
M1*(x(1)-x(3))*r12;...
M1*(x(2)-x(4))*r12];
Время работы: 8.4 с.
Производительность матлаба оказалась, на первый взгляд, на высоте, но результат несколько обескураживает:
Очевидно, что функция ode45 в Matlab требует дополнительной настройки. Для того, чтобы получить такой же результат, который дают Mathematica и фортрановская программа, пришлось несколько модифицировать матлабовский сценарий (добаленный текст выделен полужирным):
% File: orbit.m
function orbit()
global M1;
global M2;
M1=50.0; M2=10.0; V = 4.0;
tic
options=odeset('RelTol',0.0000001,'AbsTol',0.0000001,...
'InitialStep',1.0, 'MaxStep', 1000000);
[t,h]=ode45(@body2,[0,401*1000],[0,0,5,0,0,-V*M2/M1,0,V],options);
x1=h(:,1); y1=h(:,2);
x2=h(:,3); y2=h(:,4);
plot(x1,y1,x2,y2);
toc
function f=body2(t,x)
global M1;
global M2;
r12=sqrt((x(1)-x(3))^2+(x(2)-x(4))^2);
r12=1.0/r12^3;
f=[x(5);x(6);x(7);x(8);...
M2*(x(3)-x(1))*r12;...
M2*(x(4)-x(2))*r12;...
M1*(x(1)-x(3))*r12;...
M1*(x(2)-x(4))*r12];
Время расчета увеличилось до 22.6 с
Это примерно соответствует времени Mathematica и Maple. Но, по крайней мере, результат стал похож на правду:
Несколько неожиданным оказался тот факт, что функция численного решения дифференциальных уравнений требует в Matlab дополнительных настроек, в то время как в Mathematica и Maple всё корректно считается с первого раза с настройками по умолчанию. Очевидно, что при численном моделировании желательно иметь достаточно точное представление о том, что именно делает та или иная функция, и с этой точки зрения мне, пожалуй, больше подойдёт фортрановская программа, полную реализацию которой я вижу перед глазами. Да и время работы ее в три раза меньше, чем время работы матлабовского скрипта.
Ради очистки совести я запустил оба приведенных здесь матлабовских скрипта на Octave 4.0.0 (в скрипт для решения задачи двух тел в начало функции orbit() необходимо добавить строку: pkg load odepkg).
Результаты получились следующие:
Умножение матриц 5000x5000: 7.0 с.
Задача двух тел (1000 периодов): 327.3 с
Последний результат, к сожалению, делает Octave совершенно неприемлемым для моих задач. Если кто-нибудь знает, как ускорить odepkg в octave, буду благодарен за совет.
С другой стороны, при запуске первого варианта скрипта для задачи двух тел Octave сразу же выдал мне предупреждение:
warning: Option "RelTol" not set, new value 0.000001 is used
warning: called from
ode45 at line 113 column 5
orbit at line 10 column 5
warning: Option "AbsTol" not set, new value 0.000001 is used
warning: Option "InitialStep" not set, new value 40100.000000 is used
warning: Option "MaxStep" not set, new value 40100.000000 is used
С подставленными им значениями по умолчанию скрипт посчитался за 160 секунд и выдал, хотя и не вполне корректный, но, по крайней мере, похожий на правду результат:
Теперь я сведу все результаты в одну табличку:
Умножение матриц Задача двух тел
GFortran 82,3 с 7,6 с
Maple N/A ~17,0 c
Mathematica 21,9 с 16,9 с
Matlab 4,3 с 22,6 с
Octave 7,0 с 327,3 с
Я, конечно, ещё поиграюсь с матлабом до истечения триального периода, но пока склоняюсь к тому, чтобы вообще отказаться от проприетарных систем компьютерной алгебры; для решения численных задач небесной механики продолжить использовать фортран, для работы с матрицами - Octave, а для символьной математики - Maxima.