Комментарии 44
Добрые люди советуют апроксимировать дискретную функцию сплайнами, и дальше аналитически считать первую (вторую, ...) производную.
Другие добрые люди аж целую книжку на эту тему запилили. Ну или вот.
Понятно что кубическим сплайнами дальше третьей производной мы не продвинемся, но нам по условиям задачи вроде только первая нужна.
Большое спасибо за интерес к теме и за ссылки. Здесь, например, приведена формула "односторонней производной" (one-side difference), которую я и предложил в части первой для "улучшения" функции numpy.gradient на границах отрезка. Только не пойму, почему разработчики библиотек python этого не знают, или не хотят знать. А кто-то ведь им доверяет и пользуется тем, что есть.
Касательно сплайнов - при всем уважении к этому методу, требование аппроксимации функции точно проходить через узловые значения функции осложняет алгоритм. И скорее всего (я не проверял), не дает наилучшую аппроксимацию по мере 
. По другой мере, пожалуйста, только это уже другой метод. А какой лучше - попробуй определи. Скорее всего, зависит от практической задачи. В части третьей дам алгоритм ближе к практической постановке, и там реально получается кубическая интерполяция.
На прошлой неделе беседовал со своим большим другом, который лет 20 назад занимался расчетом форм для отливок из пластмассы и штампов для корпусных деталей к автомобилям. Оказывается, там нужна неперывная третья производная. Иначе отражение на поверхности автомобиля будет некрасивым.
А книжка, на которую Вы также дали ссылку, сотку баксов стоит. При всей моей любви к науке, именно на это потратить не готов. У меня, все-таки другие научные интересы.
...книжка ... сотку баксов стоит
Google Books предлагают в свободном доступе пусть и немного "усеченную", но вполне рабочую версию.
Насчет "наилучшей апроксимации" согласен, все зависит от задачи.
Спасибо огромное за продолжение, жду с нетерпением (правда :) ) третьей части.
Третья часть, наверное, будет примерно через месяц. Есть еще два-три "горячих" пирожка на другие темы. Хочу сначала их опубликовать.
Вот кстати ребята запилили некий МКЭ с участием (в том числе) B-сплайнов. Базируясь, насколько я смог понять, на разработках Micula & Micula
Понимаете, какая штука. В классическом понимании МКЭ подразумевает разделение области на небольшие части (КЭ), аппроксимацию искомой функции на каждом из этих КЭ чем-нибудь простым - линейной (простейший случай), или полиномиальной или даже тригонометрической функцией с последующей "склейкой" решений на границе соседних элементов. Если эти функции "сосредоточены" внутри одного КЭ и его соседей, то матрица СЛУ будет иметь околодиагональный вид, и СЛУ можно легко решать. Для сплайнов, помоему, получается 5 диагоналей, что уже сложнее, чем 3, но тоже терпимо. Можно развиваться в сторону усложнения этих функций. Это - один путь.
С другой стороны, если подходить к МКЭ как методу наилучших (сопряженных) аппроксимаций, то функции базиса могут быть ненулевыми хоть на всем отрезке или области. Пример. Аппроксимируем функцию на множестве полиномов. Если в качестве меры выбрать метод минимальных квадратов (а почему нет?), то Вы получите полином Лагранжа, т.е бяку. Если вы выберите 
, то аппроксимация будет вполне себе гладкой. Только матрица СЛУ получается "заполненной". Можно двигаться в этом направлении.
Можно комбинировать эти подходы и получать новые методы. Важно оченивать точность и вычислительные затраты.
На самом деле, МКЭ в одномерном случае - это из пушки по воробьям. Основное преимущество МКЭ - возможность решать 2- и 3-мерные задачи в УрЧП, причем в областях с криволинейной границей. Этим, я собственно, когда-то занимался, и сейчас постепенно вспоминаю. А вот на сплайны в 2- и 3- мерных случаях интересно было бы увидеть. Возможно, последняя ссылка об этом, там точно двумерные сетки. Посмотрю.
Был такой метод Ритца- где идея была такая- вот есть решение, давайте возьмем какой-нибудь базис функций на всей области решения, и найдем наилучшую аппроксимацию наших уравнений этим базисом. А дальше все Ваши рассуждения про L2-норму, свертки функций и минимизацию отклонений- и получалась большая плотная матрица. А потом базис из функций, заданных на всей области решения заменили на базис из финитных функций (имеющих ограниченную область поддержки, то есть, ненулевых только в каком-то небольшом объеме)- и получился МКЭ. Причем, МКЭ на симплексных сетках (одномерный МКЭ, двумерный на треугольных сетках, трехмерный на тетраэдрах) строго совпадает с конечными разностями на этих же сетках, потому что и в МКЭ и в конечных разностях базисная функция- линейная, а раз базисы одинаковые- то и аппроксимация одна и та же, хотя формулы и выглядят по разному. Но у МКЭ получилось гораздо лучше разнести операции- построение сетки- это одно, построение базиса на этой сетке- это другое, построение СЛАУ на этом базисе под эту задачу- это третье, решение СЛАУ- это четвертое. И эти операции друг с другом вообще никак не связаны, а это сильно упрощает разработку софта, реализующего вычислительную мат-физику. Выбор типа элемента никак не затрагивает построение сетки, и никак не влияет на решение СЛАУ. Физическая задача ничего не знает об элементах, она только задает функцию, которую надо минимизировать. Основное преимущество- это не возможность решать задачи в 3D- основное преимущество- это очень хорошая структурированность подхода и минимальное количество связей между отдельными операциями. и хорошо обусловленная матрица, так как элементы- конечные и матрица всегда будет разреженная. Правда, для таких матриц есть и недостатки- итерационные решения для них медленно сходятся.
С основной идеей Вашего комментария абсолютно согласен. Только МКР и МКЭ близки только в одномерном случае. МКР в 2-мерном случае можно, конечно, построить на прямоугольниках. А как на треулольниках - не очень понятно.
Кроме всего прочего, МКЭ прекрасно решает задачи в вариационной постановке, а для некоторых задач это очень полезно.
Касательно итерационных методов для 2-D и 3-D задач. Метод сопряженных градиентов дает точное (!!!) решение за N итераций, где N - число неизвестных. Я как раз в него погружаюсь :))
на треугольнике строится точно такая же разностная аппроксимация первого порядка точности- в точности как и конечноэлементная: представьте, что в одномерном случае конечная разность, деленая на шаг по координате- это на самом деле интеграл по поверхности, деленый на стягиваемый этой поверхностью объем, просто объем одномерный, а поверхность- две нульмерные стороны. это позволяет обобщить конечные разности на любое число измерений :-).
вариационность постановки задачи не связана с элементами :-). вариационные задачи сами по себе лучше решаются, чем другие постановки, потому что и требования по гладкости у них ниже, и законы сохранения в них выполняются по построению, а это для численных схем вкусный бонус. если у Вас в физической задаче энергия и импульс сохраняются- то Ваши результаты будут всегда правдоподобными. могут быть полной туфтой, но при этом все равно будут выглядеть правдоподобно, и будет сложно понять, что они тутфта :-).
про сопряженные градиенты- ну да, дает. если у Вас вычисления всегда идут с абсолютной точностью. А так-то он он не очень устойчив. я не знаю, что там у Вас за задача, но для МДТТ, гидродинамики и электродинамики от CG, BCG и их вариаций народ массово уходит на всякие CORS, TFQMR, TFQMORS и другие методы на подпространствах Крылова- у них сходимость сильно лучше, чем у CG/BCG. при такой же доказательности сходимости к точному решению. И они не требуют самосопряженности разностных операторов, что иногда (почти всегда) весьма желательно и очень ускоряет расчет. И если я правильно помню, эта сходимость строгая, если вы храните все промежуточные решения, а у Вас памяти столько нет и не будет, чтоб на хорошей сетке хранить все промежуточные решения. А без них- только медленное сползание куда-то в район решения, но без гарантий его достижения. И есть еще ряд результатов, которые намекают на то, что наилучшая сходимость среди всех итерационных методов достигается на multigrid-методах, чисто теоретически- потому что обычные итерационные завязаны на норму невязки, а невязка низкочастотных мод решения всегда очень и очень маленькая, и низкочастотные моды любыми методами ловятся крайне плохо, а мультигриды эффективно повышают частоту этой невязки на грубых сетках и ускоряют сходимость прямо феноменально. но там оператор загрубления построить не так-то вдруг, хотя есть подход algebraic multigrid, который вроде как строит наилучший оператор загрубления для любой разреженной матрицы (безотносительно сеток и вообще ее природы), но это неточно, проверять надо.
Я как раз приверженец вариационных постановок. А за подробный комментарий по-поводу итерационных методов спасибо. Буду иметь в виду.
Я когда-то давно считал 2-D задачи МДТТ (не бог весть какое достижение, естественно). Упор был на нелинейность по граничным условиям (типа "отлипания" границы от опоры), там была "фишка". Сейчас хочу воскресить наработанное для 3-D. Только вот "застрял" на обычном дифференцированиии.
Только не пойму, почему разработчики библиотек python этого не знают, или не хотят знать.
все они прекрасно знают и понимают.
The gradient is computed using second order accurate central differences in the interior points and either first or second order accurate one-sides (forward or backwards) differences at the boundaries. The returned gradient hence has the same shape as the input array.
потому что они выбрали самый быстрый и самый простой расчет градиента. и совершенно правильно сделали- для 99% задач этого хватает с головой. А кому не хватает- у тех своя голова, и они должны выбрать под свою однозначно специфическую задачу свой подходящий метод, которых миллионы, и всем не угодишь.
Сплайны - один из способов сглаживания, но далеко не единственный. МКЭ с "функциями-крышками" никакого сглаживания не даст, от кусочно-линейной функции никуда не денешься. Можно "расширить" базисные функции на на несколько соседних точек, и при кусочно-линейных базисных функциях, скорее всего, получится аналог 5-, 7- и так далее точечных схем МКР, о чем Вы также пишете. Но можно "расширенные" базисные функции сделать с непрерывными производными высших порядков, и тогда сглаживание будет так или иначе реализовано. В общем-то, сплайны вписываются в эту идеологию.
потому что совершенно разные задачи. Сплайн нужны там, где надо красиво сгладить какой-то набор данных и потом по этому сглаживанию красиво что-то изобразить. А МКЭ в силу его неявной но всегда используемой вариационной природы, обеспечивает отличное сохранение главных физических инвариантов- энергии там, массы, импульса, момента импульса и других физических инвариантов- эти инварианты гарантировано сохраняются в вариационной постановке задачи, и после дискретизации этой вариационной постановки конечными элементами тоже продолжают сохраняться. А для физики это сохранение энергии сто очков форы дает любой гладкости или красивости решений и все равно обыгрывает.
Сплайны, как и вообще-то почти любые функции, можно использовать в качестве базисных для построения МКЭ. Коллега выше дал ссылку: Вот
Да можно вообще, прости господи, ДПФ/ДВП использовать, скорость неимоверная, памяти мизер надо. Но надо данные смотреть
Расшифруйте, пожалуйста, аббревиатуру.
Быстрое дискретное преобразование Фурье /быстрое дискретное вейвлет-преобразование. Соответственно, для дифференцтрования требуется просто линейная ренормировка спектральных компонент, и обратное преобразование. Ну, для вэйвлетов чуть сложнее, но и вычислительно гибче
Я в свое время копался в этих методах, но уже подзабыл. По-моему, то и другое сильно зависят от наличия гармоник в данных (в спектре). Если гармоник нет, будет проблема со сходимостью. При этом, как я уже писал в ответах на комментарии, в 1-D случае можно много всего напридумывать. А попробуйте решить хотя бы 2-D УрЧП в области с криволинейной границей. Я на Хабре задал этот вопрос спецу по БПФ, в ответ - молчание.
Слава богу, преобразования Фурье линейны, так что сколько там измерений все равно, в том числе и в scipy.fft. структура данных неважна, но есть тонкости с фазой и нормировками.
Криволинейные границы я бы решал дополнением данных нулями и аподизацией на границах, ну если надо быстро.
Или посмотрел, не лучше ли локально в полярные/сферические координаты перейти, и кусочно по каким нибудь полиномам Лежандра или что там.
Но вообще надо для себя решить принципиальные вопросы: Вам вычислительно простой метод нужен, или есть какое-то априорное знание о характере данных, или предельно общий случай, или максимально масштабируемый, данные ограниченные или потоковые.... да много чего. У меня, разумеется, своий набор подходов, оптимальный для моих задач. В крайнем случае, надо литературу читать, наверняка Вашу узкую задачу кто-то уже пытался решить.
У меня задача - "раскочегарить" собственный мозг. Это не шутка. Я математикой занимался больше 30 лет назад.
Дело благородное. Я то числами не от хорошей жизни занимаюсь, но она , жизнь, то есть, заставляет: вот крайний раз Био-Савар в Вашем любимом МКЭ не сходился, так что дифференцировал векторные потенциалы в разрывных гранусловиях ... а казалось бы, тривиальная, школьная почти задача --всего-то посчитать магнитное поле прямоугольной катушки с током с точностью в 8 знаков.
Вейвлет хаара, обобщенный на треугольники- я делал. Если это ДВП применять на производную искомой функции, то получается, неожиданно, многосеточный метод с отличной сходимостью.
дискретное Фурье и вейвлет преобразования, наверно. только какое отношение ДПФ имеет к конечным элементам? ДПФ рулит там, где мы ОДУ во времени решаем, а МКЭ- где мы ДУ ЧП (диф.ур. в частных производных) решаем. Фурье можно иногда заюзать для решения ДУ ЧП, но в специфических случаях. В обычных практических задачах, когда коэффициенты уравнений по пространству меняются- там Фурье не дает каких-то волшебных преимуществ. А вот вариационная постановка и выполнение законов сохранения- дают.
Задача поставлена -- дифференцирование скалярного поля. ДУЧП -- это вопрос отдельный, несоразмерный поставленной задаче, и вообще бессмысленный в обобщённой формулировке. Ну, по крайней мере менее чем в 2х томах ;)
Узкая, но частная задача, где нужно дифференцировать скалярные или векторные потенциалы, известные, например, из МКЭ -- описание нетривиальных магнитных систем, где иногда требуется высокая точность.
не годится.Не очень понятно, как "обычным" способом можно продифференцировать Дельта-функцию.
Зачем "обычным", если производные известны аналитически и дискретизуются. Фокус смещён в разложение на компоненты с известными производными. Погуглите метод разложения на главные компоненты, он относительно лёгкий в вычислении, если памяти хватит, там куча частных упрощенных схем известно/можно наколхозить под задачу, типа какой-нибудь "гусеницы".
А вообще Вы уверены, что книжка Калиткина по числам не даст вам ответов на все вопросы. Раньше в мизерных выч. мощностяхиобъёмах памяти всё грамотно делали.
Если производные известны аналитически, это уже совсем другая задача, скорее к символьному дифференцированию. Но меня смущает существование (?) аналитической да еще и дискретизируемой производной Дельта-функции. Или одна обобщенная функция представляется через другие обобщенные функции? Возможно, Вы имеете в виду метод граничных элементов, когда общее решение представляется в виде комбинации фундаментальных.
Я специально несколько раз упомянул слово "разложение": да, исходя из знания типичного характера данных зачастую оптимально раскладывать в функциональные ряды, если это вычислительно легко, а разложения ортогональны или приближенно ортогональны. А аналитическая функция или обобщённая, всё равно, главное табулировать осмысленно и аккуратно.
а никак. она только в обобщенных производных дифференцируется. но от обобщенных производных один шаг до вариационных постановок :-) а от них как в черную дыру засосет в МКЭ.
Это понятно. Только, во-первых, надо еще обратно завернуть. Во-вторых, что со сходимостью - непонятно. В третьих, никто не может ответить на вопрос - как ДПФ решает задачи на областях с криволинейной границей.
Как правильно дифференцировать дискретные функции (Часть 2. Все-таки, МКЭ?)