Выбор ПО для численных методов

[chainreaction]

Осмелюсь предложить Python.
Это с одной стороны ОО язык программирования, а с другой, при подключении расширений NumPy + SciPy + MatplotLib неплохой открытый аналог Octave/Scilab.
После подключения SymPy превращается в аналог Maxima/Maple.
Имеет модуль с интерфейсом к Gnuplot, MayaVi, OpenDX и др.

Есть транслятор MatLab to Python.

http://pydev.ru/links/

Из своего горького опыта. Писал несколько месяцев в свободное время для удовольствия программу для Octave. Под самое завершение выяснилось, что без символьного дифференцирования не обойтись. Быстренько написал врапер на Python + Sympy, все работает, но пришло понимание, что все целиком изначально надо было писать на Python.

[d.dragon.n76]

octave или SciLab или MathLab удобны если студенты не знакомы с программированием на C++ или Pascal или Java ... , похож на эти средства и Fortran.

Если изучали, например, C++, то можно и изучать алгоритмы путем их непосредственного программирования, и их использовать подключая boost http://www.boost.org/doc/libs/1_40_0/libs/libraries.htm#Math . Например, решить СЛАУ методом LU разложения можно так

x=boost::math::tools::solve(A,b);

, а вычислить спецфункцию, например Бесселя, так

boost::math::cyl_bessel_k(0,x)

.
Правда возможности выполнения символьных вычислений в C++ не видел.

-- Сб ноя 14, 2009 23:18:55 --

chainreaction в сообщении #250499 писал(а):

Под самое завершение выяснилось, что без символьного дифференцирования не обойтись. Быстренько написал врапер на Python + Sympy, все работает, но пришло понимание, что все целиком изначально надо было писать на Python.

Очень интересно!

Могли бы привести пример кода программы, в котором выполняется одновременно и операция символьного дифференцирования и численная операция. Например, вычислить $\frac{d}{dx}\left(x^2\right)$ в какой-нить точке.

[Mottle]

Так подойдёт?

Используется синтаксис Python

import sympy as sp

x=sp.Symbol('x')

div=sp.diff(x**2,x,1)
sp.pprint(div)

x0=10.
print div.subs(x,x0)

 

[d.dragon.n76]

Отлично!


Помогите сделать на Питон след. программу, проверим скорость. Обратите внимание на то, как проигрывает в скорости C++.

$ time c++ slae.cc -o slae 
real   0m11.125s
user   0m4.520s
sys   0m0.412s

$ time ./slae
real   0m35.219s
user   0m33.702s
sys   0m0.012s

$ time ./slae.m 
real   0m1.678s
user   0m0.448s
sys   0m0.116s

Используется синтаксис C++

main(){
  int n = 400;
  ub::matrix<double> A(n,n);
  ub::vector<double> x(n), f(n);
  for(size_t i=0; i<n; ++i){
     for(size_t j=0; j<n; ++j){
         A(i,j) = (i==j)?10.*(0.0001*(i+j)+1):0.0001*(i+j);
     }
     f(i) = 0.01*i;
  }
  x = boost::math::tools::solve(A,f);
}
 

1. #!/usr/bin/octave -q 
2. n=400; 
3. i=[0:n-1]; 
4. [ij,ii] = meshgrid(i); 
5. A = 10.*eye(n).*(0.0001*(ij+ii)+1) .+ 0.0001*(ones(n,n).-eye(n)).*(ij.+ii); 
6. f = 0.01*i; 
7. g = A \ f'; 

[d.dragon.n76]

готовые специализированные библиотеки рулят, или шаблоны спп очень неэффективны по времени. octave проигрыват немного python наверное, из-за того, что подгружает что-то лишнее.

n=1000

$ time ./slae
real14m19.619s
user13m24.722s
sys0m0.516s

$ time ./slae.m
real0m15.288s
user0m1.032s
sys0m0.480s

$ time ./slae.py
real0m5.671s
user0m1.144s
sys0m0.396s

$ time ./slae.m
real0m8.836s
user0m1.064s
sys0m0.428s

$ time ./slae.py
real0m4.735s
user0m1.120s
sys0m0.412s

$ time ./slae.m
real0m9.214s
user0m1.048s
sys0m0.436s

Используется синтаксис Python

#!/usr/bin/python
import numpy as np
n=1000
i=np.linspace(0,n-1,n)
[ij,ii]=np.meshgrid(i,i)
A=10*np.eye(n)*(0.0001*(ij+ii)+1)+0.0001*(np.ones((n,n))-np.eye(n))*(ij+ii)
f = 0.01*i
x=np.linalg.solve(A,f)

 

[nestoklon]

d.dragon.n76 в сообщении #262077 писал(а):

Например, решить СЛАУ методом LU разложения можно так
Код:
x=boost::math::tools::solve(A,b);

http://alxr.usatlas.bnl.gov/lxr/source/external/Boost_1.44.0_python2.6/boost-1_44/boost/math/tools/solve.hpp писал(а):

033 //
034 // Caution: this uses undocumented, and untested ublas code,
035 // however short of writing our own LU-decompostion code
036 // it's the only game in town.
037 //

Причём настолько undocumented, что на сайте буста этой функции нет.

Не надо пользоваться линейкой из буста. Она там плохая. И из GSL линейкой тоже пользоваться не надо. Надо как все нормальные питоны делают вызывать Lapack.

[d.dragon.n76]

nestoklon в сообщении #427008 писал(а):

Причём настолько undocumented, что на сайте буста этой функции нет.

boost тоже не стандарт относительно STL, например .

solve это всего лишь интерфейс к LU-разложению. Можно использовать /usr/include/boost/numeric/ublas/lu.hpp, там сказано "// LU factorizations in the spirit of LAPACK and Golub & van Loan", и задокументированы эти функции. Я кроме чтения док. просматриваю /usr/include/... и заметил полезные модули.

Встречали ли вы интерфес на C++ к LAPACK? По моему, желательно пользоваться чем-то общеизвестным, типа STL или boost. На boost.org вижу только интерфейс к BLAS --- это uBLAS (на ЧМ, кстати, вполне можно использовать при программировании ЧМ использующих матричные операции), что-то типа LAPACK для С++ почему-то не написали. С API, да и C считаю тяжеловатыми при решении задач, в которых не нужно сосредотачиваться на программировании классического ЧМ.

nestoklon в сообщении #427008 писал(а):

Не надо пользоваться линейкой из буста. Она там плохая.

А вот здесь обоснуйте, хочу разобраться. Например, почему BLAS из boost хуже таких же функций из Питон?

[nestoklon]

d.dragon.n76 в сообщении #427230 писал(а):

solve это всего лишь интерфейс к LU-разложению.

Простите, Вы читаете что Вам пишут? Это нетестированный и недокументированный код.

d.dragon.n76 в сообщении #427230 писал(а):

Встречали ли вы интерфес на C++ к LAPACK?

Его забросили так и не доделав. Lapack++ называется. Предполагаю, никто не пользовался.

d.dragon.n76 в сообщении #427230 писал(а):

С API, да и C считаю тяжеловатыми при решении задач, в которых не нужно сосредотачиваться на программировании классического ЧМ.

Клгда-то давно долго искал. Лучшее что нашёл эта библиотека. Но потом привык к мысли, что проще тупо вызывать Lapack. Благо, он настолько стандартный, что даже названия функций стандартизированы между различными библиотеками лет много как.

d.dragon.n76 в сообщении #427230 писал(а):

Например, почему BLAS из boost хуже таких же функций

Не знаю. Это опытный факт. Программисты писали?.. (Это шутка конечно, BLAS там вряд ли отличается, вот только Вы в курсе что LU разложения в BLAS нет?)

[d.dragon.n76]

nestoklon в сообщении #427249 писал(а):

d.dragon.n76 в сообщении #427230 писал(а):

solve это всего лишь интерфейс к LU-разложению.

Простите, Вы читаете что Вам пишут?

Читаю, конечно же. Шаблон solve не тестированный и не задокументированный, согласен. Я говорил о том, что всегда есть возможность воспользоваться документированными возможностями. И отмечаю то, что библиотеки boost это кандидаты на стандартизацию ("не до STL" ), и с тех позиций что кто-то не протестировал и не описал, чем то кроме STL в cpp пользоваться опасно.


В документации LU-разложение я найти его не могу, хотя в /usr/include/boost/numeric/ublas/lu.hpp есть шаблон:

Используется синтаксис C++

// LU factorization with partial pivoting
 template<class M, class PM>
  typename M::size_type lu_factorize (M &m, PM &pm) 
 

[nestoklon]

d.dragon.n76 в сообщении #427256 писал(а):

и с тех позиций что кто-то не протестировал и не описал, чем то кроме STL в cpp пользоваться опасно

Много лет назад мне рассказали сказочную историю про стандартные библиотеки. Один знакомый программист долго искал в своей программе ошибку. Очень долго искал. В конце концов начал просматривать что происходит не только в его коде, но и в стандартных библиотеках (а пользовался он стандартной библиотекой идущей в комплекте с VS). В итоге нашёл в шаблоне sort код приблизительно следующего содержания: "если длина вектора больше чем <некоторое достаточно большое чисто>, возвращаем без сортировки".
Буквально на днях при попытке понять где у меня утекает память, нашёл в описании библиотеки которой пользуюсь, сказочную рекомендацию -- не пользоваться на 64 битном линуксе библиотекой glibc, потому то в ней есть баг с выделением памяти, а подменять её аналогичной библиотекой от гугла (которая, как утверждается, к тому же быстрее).
Это я к чему? Не надо рассчитывать на то, что если что-то стандартно, то оно обязательно работает правильно и оптимально. boost -- это своего рода обкатка новых фич. Они конечно заботятся о скорости, эффективности и прочем, но у них не всегда получается. О том, чтобы оно не ломалось -- они заботятся по-настоящему (причём неизвестно, насколько этой цели добиваются), но не о скорости точно. Тем более что большие куски его устаревают когда вдохновлённые этой обкаткой авторы стандарта решают добавить поддержку некоторых из этих фич в язык. А Вы берёте недокументированную фичу и удивляетесь, а чего это она хреново работает? Потому и недокументированная.

[d.dragon.n76]

nestoklon писал(а):

...

Что вы предлагаете использовать студентам в курсовых и дипломных, например, по математическому моделированию? На программировании студенты изучали cpp (или pascal). Преподаватели численных методов, и других близких дисциплин, имеют возможность в выборе библиотек для изучаемых языков или, вообще, других языков (скриптовых, специализированных) и заинтересованы в том, чтобы студент справился с под.задачами из научных работ, типа численно решить СЛАУ, нарисовать график, численно посчитать интеграл и другое.

[nestoklon]

d.dragon.n76 в сообщении #427370 писал(а):

Что вы предлагаете использовать студентам в курсовых и дипломных, например, по математическому моделированию?

Питон.

[myhand]

nestoklon в сообщении #427008 писал(а):

И из GSL линейкой тоже пользоваться не надо.

А с GSL-то что тут не так?

[GrayDog]

Я использую SciLab. Свободный, кроссплатформенный и достаточно мощный пакет.

[Snegovik]

Современные версии Фортрана тоже вполне актуальны.