Хабр Курсы для всех
РЕКЛАМА
Практикум, Хекслет, SkyPro, авторские курсы — собрали всех и попросили скидки. Осталось выбрать!
Хабр доступен 24/7 благодаря поддержке друзей
Комментарии 24
Лет так двадцать-тридцать назад для решения какого-нибудь нетривиального ОДУ нужно было знать кучу методик из вычматов, разных аппроксимаций и уточнений, сейчас же все в одну строку одним методом делается. Человечество далеко шагнуло вперед в этой области, спору нет.
для решения какого-нибудь нетривиального ОДУ нужно было знать кучу методик из вычматов, разных аппроксимаций и уточнений
Оно и сейчас бывает нужно, смотря какая задача — у меня недавно целая история была с жесткой системой, да и не уверен что распутал. А в целом да — типовые задачи решаются в пару строк
Своя, отечественная система это замечательно.
Но, вот например, только Python-овские пакеты sympy, scipy (с numpy и matplotlib + внешние, тот же GSL) давно делают, то же самое, что SMath Studio. Даже намного больше.
Если нужен интерфейс, то опять Python-овский sage.
И главное, они давно отлажены, они бесплатные и в них и для них работают многие.
Я просто знаю, что разработка нормальной (универсальной) для работы системы компьютерной алгебры занимает примерно 100 человеко-лет.
Ниже примеры, которые несколько лет назад я писал для студентов.
Но, вот например, только Python-овские пакеты sympy, scipy (с numpy и matplotlib + внешние, тот же GSL) давно делают, то же самое, что SMath Studio. Даже намного больше.
Если нужен интерфейс, то опять Python-овский sage.
И главное, они давно отлажены, они бесплатные и в них и для них работают многие.
Я просто знаю, что разработка нормальной (универсальной) для работы системы компьютерной алгебры занимает примерно 100 человеко-лет.
Ниже примеры, которые несколько лет назад я писал для студентов.
Уравнение Ван-дер-Поля (представление Льенара)
# -*- coding: utf-8 -*-
"""
Уравнение Ван-дер-Поля (представление Льенара)
\begin{gather*}
y^{\prime}_1 = - a (\frac{y_1^3}{3} - y_1) + a y_2, \\
y^{\prime}_2 = - y_1,
\end{gather*}
Здесь a = 103 и a = 106, y1 = 2, y2= 0.
t_k = 20.
"""
from sympy import *
from scipy.integrate import ode
from pylab import *
y1, y2 = Symbol('y[0]'), Symbol('y[1]')
a = Symbol('a')
f1 = -a*(y1**3/3 - y1) + a*y2
f2 = -y1
Y = Matrix([f1, f2])
X = Matrix([y1, y2])
print 'f1 =', f1
print 'f2 =', f2
print 'jacobian = ', Y.jacobian(X)
def f(t, y, a):
return [a*y[1] - a*(-y[0] + y[0]**3/3),
-y[0]]
def jac(t, y, a):
return [[-a*(1 - y[0]**2), a],
[-1, 0]]
y01, y02 = 2.0, 0.
a = 103.
r = ode(f, jac).set_integrator('vode', method='bdf', with_jacobian=True)
r.set_initial_value((y01, y02)).set_f_params(a).set_jac_params(a)
t1, dt, t, y1, y2 = 20, 0.1, [0.], [y01], [y02]
while r.successful() and r.t < t1:
r.integrate(r.t + dt)
t.append(r.t)
y1.append(r.y[0])
y2.append(r.y[1])
figure(figsize=(13, 18))
splt1 = subplot(3, 1, 1)
title("$y^{\prime}_1 = - a (\\frac{y_1^3}{3} - y_1) + a y_2,\, y^{\prime}_2 = - y_1,\, a = %.2f,\, y1(0) = %.2f, \,y2(0)=%.2f$" % (a, y01, y02))
ylabel(r'$y_1$', {'fontsize': 18})
xlabel(r'$t$', {'fontsize': 18})
plot(t, y1, 'o-')
grid(True)
splt2 = subplot(3, 1, 2)
ylabel(r'$y_2$', {'fontsize': 18})
xlabel(r'$t$', {'fontsize': 18})
grid(True)
plot(t, y2, 'o-')
splt3 = subplot(3, 1, 3)
ylabel(r'$y_1$', {'fontsize': 18})
xlabel(r'$y_2$', {'fontsize': 18})
grid(True)
plot(y2, y1, 'o-')
savefig('./Уравнение Ван-дер-Поля.pdf')
show()
Результат
f1 = a*y[1] — a*(y[0]**3/3 — y[0])
f2 = -y[0]
jacobian = Matrix([[-a*(y[0]**2 — 1), a], [-1, 0]])
f2 = -y[0]
jacobian = Matrix([[-a*(y[0]**2 — 1), a], [-1, 0]])
Уравнение Бонгоффера - Ван-дер-Поля
# -*- coding: utf-8 -*-
"""
Уравнение Бонгоффера - Ван-дер-Поля
\begin{gather*}
y^{\prime}_1 = a(- (\frac{y_1^3}{3} - y_1 ) + y_2), \\
y^{\prime}_2 = - y_1 - by_2 + c
\end{gather*}
Здесь a = 103 и a = 106, y1 = 2, y2= 0.
Уравнение описывает протекание тока через клеточную мембрану.
Постоянная компонента тока c в безразмерной записи системы такова,
что 0 < c < 1, b > 0. t_k = 20.
"""
from sympy import *
from scipy.integrate import ode
from pylab import *
y1, y2 = Symbol('y[0]'), Symbol('y[1]')
a, b, c = Symbol('a'), Symbol('b'), Symbol('c')
f1 = a*(-(y1**3/3 - y1) + y2)
f2 = -y1 - b*y2 + c
Y = Matrix([f1, f2])
X = Matrix([y1, y2])
print 'f1 =', f1
print 'f2 =', f2
print 'jacobian = ', Y.jacobian(X)
def f(t, y, a, b, c):
return [a*(y[0] + y[1] - y[0]**3/3),
c - y[0] - b*y[1]]
def jac(t, y, a, b, c):
return [[a*(1 - y[0]**2), a],
[ -1, -b]]
y01, y02 = 2., 0.
a, b, c = 103., 2.1, 0.5
r = ode(f, jac).set_integrator('vode', method='bdf', with_jacobian=True)
r.set_initial_value((y01, y02)).set_f_params(a, b, c).set_jac_params(a, b, c)
t1, dt, t, y1, y2 = 20, 0.5, [0.], [y01], [y02]
while r.successful() and r.t < t1:
r.integrate(r.t + dt)
t.append(r.t)
y1.append(r.y[0])
y2.append(r.y[1])
figure(figsize=(13, 18))
splt1 = subplot(3, 1, 1)
title("$y^{\prime}_1 = a (- (\\frac{y_1^3}{3} - y_1 ) + y_2),\, y^{\prime}_2 = - y_1 - b y_2 + c,\, a = %.2f,\, b = %.2f,\, c = %.2f,\, y1(0) = %.2f, \,y2(0)=%.2f$" % (a, b, c, y01, y02))
ylabel(r'$y_1$', {'fontsize': 18})
xlabel(r'$t$', {'fontsize': 18})
plot(t, y1, 'o-')
grid(True)
splt2 = subplot(3, 1, 2)
ylabel(r'$y_2$', {'fontsize': 18})
xlabel(r'$t$', {'fontsize': 18})
grid(True)
plot(t, y2, 'o-')
splt3 = subplot(3, 1, 3)
ylabel(r'$y_2$', {'fontsize': 18})
xlabel(r'$y_1$', {'fontsize': 18})
grid(True)
plot(y1, y2, 'o-')
savefig('./Уравнение Бонгоффера - Ван-дер-Поля.pdf')
show()
Результат
f1 = a*(-y[0]**3/3 + y[0] + y[1])
f2 = -b*y[1] + c — y[0]
jacobian = Matrix([[a*(-y[0]**2 + 1), a], [-1, -b]])
f2 = -b*y[1] + c — y[0]
jacobian = Matrix([[a*(-y[0]**2 + 1), a], [-1, -b]])
С полной поддержкой единиц измерения (локализованных)? С интерфейсом аля Mathcad?
SMath Studio тоже бесплатна и расширяемая.
SMath Studio тоже бесплатна и расширяемая.
Единицы измерения в sympy есть, правда я пользовался ими только в достаточно простых случаях. Под локализацией вы подразумеваете перевод названий на русский? Не думаю, что есть смысл — единицы измерения и так понятны, но вряд ли это сложно добавить.
А из Mathcadовского интерфейса что именно вам нужно? Всякие панельки для ввода символов мышкой намного медленнее, чем ввод с клавиатуры, а объединение в одном документе кода, текста (+latex), графиков и всего такого — это отличнейшим образом сделано в IPython notebook.
В общем, если не хватает чего-то небольшого вроде единиц измерения, то намного продуктивнее было бы дописать их к sympy, а не разрабатывать всё с нуля. А преимущества Python очевидны — это всё-таки полноценный язык программирования, можно при необходимости добавить что угодно, да и библиотек полно.
А из Mathcadовского интерфейса что именно вам нужно? Всякие панельки для ввода символов мышкой намного медленнее, чем ввод с клавиатуры, а объединение в одном документе кода, текста (+latex), графиков и всего такого — это отличнейшим образом сделано в IPython notebook.
В общем, если не хватает чего-то небольшого вроде единиц измерения, то намного продуктивнее было бы дописать их к sympy, а не разрабатывать всё с нуля. А преимущества Python очевидны — это всё-таки полноценный язык программирования, можно при необходимости добавить что угодно, да и библиотек полно.
Под локализацией я понимаю полную поддержку единиц измерения на всех языках сразу (как сделано в SMath Studio), если разночтение имеет место. Если я буду делать расчёты в школе, то мне что, ещё и английский изучать, чтобы единицы измерения писать? А зачем? Если в учебнике всё по-русски, то и оформлять я тоже хочу по-русски. А если расчётный документ нужно приложить официально с использованием размерностей, то зачем мне там импортные эквиваленты?
Я не против Python'а, даже наоборот, я включил поддержку скриптов на c#, vb.net и python, прямо в документе SMath Studio, написав соотв-щее дополнение к программе. Но мне нравится математический вид документа сразу, а не после компиляции, мне нравится возможность «разбросать» формулы по листу как я хочу, не «ковыряясь» в разметке. Ничего не имею против LaTeX'а, но он статичен.
В SMath Studio среда предлагает доступные функции (единицы измерения) после первого введённого символа, т.е. мышкой кликать нужно не так уж и часто.
Я не против Python'а, даже наоборот, я включил поддержку скриптов на c#, vb.net и python, прямо в документе SMath Studio, написав соотв-щее дополнение к программе. Но мне нравится математический вид документа сразу, а не после компиляции, мне нравится возможность «разбросать» формулы по листу как я хочу, не «ковыряясь» в разметке. Ничего не имею против LaTeX'а, но он статичен.
В SMath Studio среда предлагает доступные функции (единицы измерения) после первого введённого символа, т.е. мышкой кликать нужно не так уж и часто.
Про единицы измерения ещё раз напишу — локализацию их, я думаю, добавить в sympy несложно. Но всё-таки, разве есть кто не поймёт, что m это метр, s — секунда и другие общеупотребительные величины? Сложные и редко используемые и по русской аббревиатуре бывает непонятно.
Насчёт «разбросать по листу» в 2D — именно такого из коробки в питоне нет, но вот организовать линейную структуру с секциями легко и удобно можно в IPython Notebook — вы попробуйте) Latex я упомянул только в плане формул, то что их можно вставлять прямо в Notebook.
Я это всё к чему — сам пробовал когда-то SMath Studio, но не показалось удобно для чего-то кроме простейших примеров. В обычной школе показывать я бы тоже выбрал её, наверное, но вот дальше — только полноценные языки программирования, где точно можно сделать всё что нужно.
Насчёт «разбросать по листу» в 2D — именно такого из коробки в питоне нет, но вот организовать линейную структуру с секциями легко и удобно можно в IPython Notebook — вы попробуйте) Latex я упомянул только в плане формул, то что их можно вставлять прямо в Notebook.
Я это всё к чему — сам пробовал когда-то SMath Studio, но не показалось удобно для чего-то кроме простейших примеров. В обычной школе показывать я бы тоже выбрал её, наверное, но вот дальше — только полноценные языки программирования, где точно можно сделать всё что нужно.
SMath Studio на самом деле ещё далека до реального использования где-то кроме школы. Делает практически один человек.
У нас в институте преподаватель один не принимал расчёты без единиц измерения, считая их безграмотными. Пусть цифры правильные, но пока не приведёшь всё к размерному виду — работу не зачтёт.
Думаю, что и официальные органы, где нужно выполнять расчёты по СНиПам, ГОСТам и прочим нормативным документам, не оценят импортные размерности. Посчитать вы можете и с ними, но вот в отчёте всё нужно будет привести к отечественным эквивалентам. Я как-то делал такой расчёт.
У нас в институте преподаватель один не принимал расчёты без единиц измерения, считая их безграмотными. Пусть цифры правильные, но пока не приведёшь всё к размерному виду — работу не зачтёт.
Думаю, что и официальные органы, где нужно выполнять расчёты по СНиПам, ГОСТам и прочим нормативным документам, не оценят импортные размерности. Посчитать вы можете и с ними, но вот в отчёте всё нужно будет привести к отечественным эквивалентам. Я как-то делал такой расчёт.
Вот простой пример.
Я бы добавил в этот список Моделику.
Меня больше интересовали законченные решения в виде готовых библиотек, которые можно было бы применять в своих программах на ЯВУ. Математические пакеты хороши для прототипирования, математического моделирования, в общем тогда, когда есть время порешать, подумать, а вот когда всё рассчитано и проверено, хочется быстродействия и нужно переписать всё действо в другом виде, с использованием окошек, кнопочек и пр., то встаёт вопрос где брать эквиваленты математических операций из пакетов.
По этой причине мат. пакеты обзавелись встроенными возможностями по созданию интерфейсов или по переводу программ в компилируемый вид. Но такая «поддержка» требует иметь довесок в несколько сотен мегабайт, например, у Matlab последних версий.
Вот в общем и весь смысл, мне хотелось иметь небольшие библиотечки, которые я бы мог безболезненно прицепить к программе на c#, с++ или delphi. Решатели ОДУ — вещь специфическая, а жёстких ОДУ тем более.
По этой причине мат. пакеты обзавелись встроенными возможностями по созданию интерфейсов или по переводу программ в компилируемый вид. Но такая «поддержка» требует иметь довесок в несколько сотен мегабайт, например, у Matlab последних версий.
Вот в общем и весь смысл, мне хотелось иметь небольшие библиотечки, которые я бы мог безболезненно прицепить к программе на c#, с++ или delphi. Решатели ОДУ — вещь специфическая, а жёстких ОДУ тем более.
Для R есть замечательный пакет deSolve, в котором находится большое количество различных численных методов для решения ОДУ. Мой любимый — lsoda, он умеет автоматически переключаться между жёсткими и нежёсткими системами. По сути это обёртка над FORTRAN-библиотекой за авторством Линды Петзольд и Алана Хиндмарша из Ливерморской национальной лаборатории. Библиотека начала разрабатываться ещё в 70-ые, текущая версия базируется на имплементации 2003-года из библиотеки NetLib.
Рекомендую потестировать библиотеку SADEL: pa10.ru/
Это интересно, спасибо. Тестирование там уже проведено, но я бы тоже не отказался попробовать решатель. Что-то пока зарегистрироваться не могу, чтобы скачать библиотеки. Как протестирую, включу в обзор.
А вы метод Адамса не используете? Мне он показался более быстрым и устойчивым. Кроме этого там проще осуществить контроль шага
Разве метод адамса-башфорта или адамса-мултона предусматривают контроль шага? Мне казалось, что данные методы, вернее их аппроксимирующие полиномы, построены на принципе постоянства шага?
Вот самописный maple-код для генерации формул методов Адамса любого порядка. Порядок задается параметром p
В плане контроля локальной погрешности мог бы порекомендовать метод Рунге-Кутты с модификацией по Фельдбергу. Может работать и для умеренно жестких задач
Вот самописный maple-код для генерации формул методов Адамса любого порядка. Порядок задается параметром p
В плане контроля локальной погрешности мог бы порекомендовать метод Рунге-Кутты с модификацией по Фельдбергу. Может работать и для умеренно жестких задач
В некоторых задачах вместо жестких решателей можно использовать и другие, используя более мелкую сетку. Просто считать они будут дольше, если решение сойдётся. Главное, чтобы было из чего выбирать, а там уже практика покажет какой из методов даёт устойчивый результат.
Есть ещё boost::odeint, вдруг кому эта информация окажется полезной.
Добавил три ссылки: boost::odeint, SADEL и решатель Лимонова А. Г.
Автор последнего решателя обратился лично. Может кто-нибудь пригласит его на хабр? У меня похоже нет такой возможности. Было бы интересно услышать комментарии от автора диссертации на эту тему (почта: a.limonov на gmail).
Будет время, протестирую все остальные библиотеки. Пока не знаю как грамотно делать количественные сравнения жестких решателей, может быть осилю примеры из диссертации, тогда можно сделать сравнительную статью.
Автор последнего решателя обратился лично. Может кто-нибудь пригласит его на хабр? У меня похоже нет такой возможности. Было бы интересно услышать комментарии от автора диссертации на эту тему (почта: a.limonov на gmail).
Будет время, протестирую все остальные библиотеки. Пока не знаю как грамотно делать количественные сравнения жестких решателей, может быть осилю примеры из диссертации, тогда можно сделать сравнительную статью.
Только прочитал этот пост. Спасибо, буду его иметь в виду при решении своих задач. Возникли вопросы:
1. Как известно, для численного решения жестких ОДУ используются неявные разностные схемы. Их недостатком является то, что в общем случае на каждом шаге интегрирования приходится решать нелинейную систему алгебраических уравнений, а она может иметь более одного решения (есть примеры). Как с этим делом справляются пакеты?
2. В литературе по современным численным методам для устойчивости неявных схем на величину шага накладываются ограничения (связь с константой Липшица). Как я понимаю, это здесь не проверяется?
1. Как известно, для численного решения жестких ОДУ используются неявные разностные схемы. Их недостатком является то, что в общем случае на каждом шаге интегрирования приходится решать нелинейную систему алгебраических уравнений, а она может иметь более одного решения (есть примеры). Как с этим делом справляются пакеты?
2. В литературе по современным численным методам для устойчивости неявных схем на величину шага накладываются ограничения (связь с константой Липшица). Как я понимаю, это здесь не проверяется?
Зарегистрируйтесь на Хабре, чтобы оставить комментарий
Обзор доступных библиотек для численного решения жёстких ОДУ