apximhd | Варианты выбора

Подходит к концу моя студенческая лицензия на 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.

Flat | Top-Level Comments Only

From:

evgeniirudnyi.livejournal.com

Мне в свое время (9 лет назад, когда она у меня была) нравилась Математика. В ней было хорошее совмещение символьных и численных вычислений. Кстати, в Математике должны быть функции с использованием LAPACK. По всей видимости это происходит не автоматически, по всей видимости надо преобразовать матрицу в какое-нибудь другое представление.

Из бесплатных программ могу рекомендовать numpy и scipy в Python. У меня кстати есть тест с перемножением матриц Matrix Multiplication на matrixprogramming com.

From:

al-pas.livejournal.com

Мне тоже нравится Mathematica, и я её достаточно давно использую. Но встал вопрос о продлении лицензии за собственный счёт, и я задался вопросом: если я буду отдавать свои кровные деньги, то продолжить ли использовать Mathematica или имеет смысл отдать деньги за что-нибудь другое. Нужна символьная математика, числодробилка и презентационные возможности в одном флаконе. Вот как-то пока альтернативы тому, что я использую сейчас, то есть, связке Mathematica + Fortran не видится, но Mathematica + Fortran это два флакона.

From:

evgeniirudnyi.livejournal.com

В Mathematica можно вызывать внешние подпрограммы через Mathlink. У меня есть пример для вызова quadpack. Также webMathematica была неплоха. Далее можно писать книги с красивыми формулами.

Однако вопрос денег как всегда не простой. Сложно сказать.

From:

Варианты выбора

Navigation

Page Summary

Style Credit

Expand Cut Tags