Часто задаваемые вопросы

Во-первых, сохраняйте спокойствие. Вы могли получить одно из следующих предупреждений:

Все эти сообщения указывают на одно и то же: по той или иной причине решение ОДУ расходится к бесконечности. Из-за расхождения к бесконечности dt интегратора уменьшается (с целью контроля скорости и погрешности) и поэтому достигает минимального значения dt, достигает dt=NaN, либо значение в ОДУ достигает Inf. В зависимости от того, что происходит первым, выдается соответствующее предупреждение.

Как быть в этом случае? В 99,99 % случаев в результате отладки выяснялось, что ошибка кроется в модели пользователя! Пропущенный знак «минус», неверный термин и т. д. Нужно обращать внимание на это и многое другое. В некоторых ОДУ увеличение значения параметра может привести к бифуркации, из-за чего решение расходится. Если в уравнении u'=a*u переменная a отрицательная, результат прекрасно сходится к нулю, но если переменная a положительная, решение быстро расходится к бесконечности! То есть нужно внимательно следить за правильностью знаков перед параметрами!

Обратите внимание: если эти предупреждения появляются в процессе оценки параметров, то, скорее всего, проблема кроется в нем. Просто проверьте SciMLBase.successful_retcode(sol.retcode) и задайте стоимость Inf. Большинство оптимизаторов не станут выполнять шаги в таких режимах параметров!

Следует проверить и еще ряд моментов. Зачастую устойчивость решения ОДУ повышается с уменьшением погрешности, поэтому можно попробовать уменьшить значения abstol и reltol. Если ваша модель представляет собой дифференциально-алгебраическое уравнение (ДАУ) с высоким индексом (например, больше 1), это может повлиять на численное решение. В таком случае можно воспользоваться инструментом снижения индекса ModelingToolkit.jl, чтобы повысить устойчивость численного решения. Кроме того, если вы используете безматричный решатель (такой как GMRES) для решения большого ОДУ или ДАУ с высокой степенью жесткости, убедитесь в том, что погрешность GMRES настроена правильно и применяется соответствующий предобусловливатель. Наконец, попробуйте другие решатели. У них у всех разная устойчивость, поэтому попробуйте использовать Tsit5(), Vern7(), QNDF(), Rodas5(), TRBDF2(), KenCarp4(), Sundials.CVODE_BDF() и т. д. и посмотрите, что работает лучше.

Если ничего не помогло, еще раз проверьте, действительно ли ваше ОДУ ведет себя так, как вы от него ожидаете. Это одна из самых распространенных причин проблемы: ваша интуиция может вас подводить. Например, u' = -sqrt(u) при u(0)=1 не может достичь нуля, потому что производная стремится к нулю, так? Нет, не так! Нуль будет достигнут через конечное время, после чего решение станет неопределенным, то есть действительного решения не будет. u' = u^2 - 100u «обычно» стремится к нулю, но при u(0)>10 стремится к бесконечности. Постройте график расходящегося решения и посмотрите, верна ли асимптотика: при увеличении u[i] становится ли производная u'[i] в уравнении положительной и растет ли ее значение? Проблема может быть в этом!

Предположим, вы уверены, что у вас ошибки нет, и хотите сообщить о проблеме. Для этого сначала нужно доказать, что она связана именно с решателем. Если проблема в решателе, то она не должна возникать во всех решателях. Считаете, что проблема в решателях Julia? На этот случай DifferentialEquations.jl предлагает прямые неизмененные оболочки почти для всех ранее созданных решателей, поэтому если вы думаете, что проблема в Julia, попробуйте пропустить ваше ОДУ через следующие оболочки:

 — LSODA.jl, оболочка для классического алгоритма lsoda;

И этим список не ограничивается. Для тестирования достаточно изменить solve(prob,Tsit5()) на solve(prob,lsoda()), так что попробуйте. Если вы перевели код с другого языка, например Python или MATLAB, используйте прямую оболочку, чтобы проверить, не изменились ли шаги. Если изменились, значит изменилось и ОДУ, потому что оно использует прямой вызов решателей из этих пакетов!

Если ОДУ расходится к бесконечности во всех мыслимых решателях, то проблема, скорее всего, не в них. Здесь будет уместен такой мем:

Не будьте как Патрик. Если вы перепробовали все рекомендации, но с решением ОДУ по-прежнему есть проблема и установить ее причину не получается, смело задавайте вопрос на форуме обсуждения Julia, чтобы получить помощь в диагностике. Если вы все же нашли проблему в решателе, сообщите о проблеме в репозитории GitHub.

Вы можете получить такое сообщение:

Имейте в виду: вполне вероятно, что причина просто в слишком длительном временном интервале. Если при проверке sol.t в возвращенном объекте видно, что величина шага приемлемая, просто передайте параметр maxiters=... в решатель со значением больше значения по умолчанию Int(1e5).

Но если значение maxiters уже достаточно велико, то, скорее всего, проблема в жесткости вашей модели. Жесткое ОДУ требует очень небольших временных шагов при использовании многих явных решателей, таких как Tsit5(), Vern7() и других, поэтому эти методы не подходят для такой задачи. Лучше выбрать другой метод, например Rodas5(), Rosenbrock23(), TRBDF2(), KenCarp4() или QNDF().

Для этого есть множество инструментов! Однако сперва давайте остановимся на таком вопросе: когда вы говорите «должно», что вы под этим подразумеваете? Если вы имеете в виду «я могу математически доказать, что данное ОДУ при таких значениях и начальных условиях должно иметь решение, которое будет положительным всегда», тогда все верно. Если же вы имеете в виду «это модель биохимических реакций, и концентрация всегда должна быть положительной», то сначала спросите у себя, правильно ли вы сформулировали модель, так чтобы решение всегда было положительным?

Описанный ниже набор инструментов предназначен для обеспечения положительности моделей на основе ОДУ, истинное решение которых согласно математическому обоснованию должно быть положительным. При обнаружении модели, которая дает отрицательное решение, они постараются найти положительное, но при этом верное решение, что невозможно, и поэтому они просто выдадут ошибку. Здесь все не так просто, как кажется. Решение u'=-sqrt(u) не обязательно остается положительным, несмотря на то, что производная стремится к нулю при стремлении u к нулю (если вам интересно, можете проверить аналитическое решение). Точно так же при анализе нелинейных моделей может выявляться самое разное поведение. Распространенный случай непредвиденного возникновения отрицательности — функции Хилла в моделях биологических систем: то, что производная стремится к нулю, не означает, что это происходит достаточно быстро, чтобы решение в целом оставалось положительным!

Учитывая, давайте посмотрим, какие есть варианты.

Простейший прием — изменение погрешности решателя. Немного уменьшите значение abstol (и, возможно, reltol). Это может помочь сократить погрешность и, таким образом, сделать решение положительным. В случае с некоторыми более сложными уравнениями может помочь выбор жесткого решателя ОДУ, такого как Rosenbrock23() QNDF или TRBDF2().

Если это не помогло, пора прибегнуть к тяжелой артиллерии. Один из вариантов — параметр isoutofdomain, позволяющий определить логическую функцию, которая будет отклонять шаги при невыполнении условий. Например, если указать isoutofdomain = (u,p,t)->any(x->x<0,u), решатель будет отклонять все шаги, приводящие к отрицательному значению переменной u. Используя любой решатель Julia с этим параметром, невозможно получить отрицательный результат! Однако вы можете получить такое предупреждение:

либо

Это означает, что обеспечить положительность невозможно. Шаг с отрицательным результатом отклоняется, значение dt уменьшается, выполняется следующий шаг, отклоняется, снова уменьшается значение, и так повторяется до тех пор, пока dt не достигнет dtmin или не будет достигнуто значение maxiters. Это означает, что даже при попытке решить задачу с наиболее точной бесконечно малой величиной dt решение все равно становится отрицательным. Уверены ли вы, что истинное решение должно быть положительным? Попробуйте поискать ошибки в уравнениях, например пропущенный знак «минус».

Если такой способ работает, но немного медленно, в библиотеке обратных вызовов есть обратные вызовы для проверки области значений, которые делают то же самое, но с большей производительностью. Вместо повторения множества шагов из-за отклонения производится обратная интерполяция, чтобы хоть небольшой, но шаг вперед был сделан. Однако устойчивость может немного снизиться, поэтому применимость данного способа зависит от уравнения, и опять-таки решение действительно должно быть положительным. Если решение становится отрицательным, обратная интерполяция будет выполняться до возможного предела, после чего произойдет ошибка, связанная с dtmin.

Наконец, имейте в виду, что правильность решения, выдаваемого решателями ОДУ, не превышает погрешность. Поэтому если решение должно быть положительным, но abstol=1e-12, в результате можно получить u[i]=-1e-12. И это вполне ожидаемое поведение численных решателей ОДУ, ведь они так или иначе выполняют поставленную задачу. Если для вашего приложения это серьезная проблема, вам может потребоваться предусмотреть такое поведение в своей модели, например изменить sqrt(u[i]) на sqrt(max(0,u[i])). Можно также преобразовать значения, например использовать выражение u^2 или exp(u) вместо u, что математически может давать только положительное значение. Обратите внимание на такие инструменты, как ModelingToolkit.jl, для автоматического преобразования уравнений.

Да. Все библиотеки *DiffEq.jl (OrdinaryDiffEq.jl, StochasticDiffEq.jl и DelayDiffEq.jl) являются универсальными по отношению к типам массивов и чисел. Это означает, что они принимают тип массива, заданный в реализации. Выполняемые на месте алгоритмы на внутреннем уровне используют механизм трансляции Julia (за некоторыми исключениями из-за имеющейся пока в Julia ошибки, см. эту проблему) и функцию Julia mul! для перемножения матриц на месте. Алгоритмы, выполняемые не на месте, используют стандартные арифметические функции. И в том, и в другом случае дополнительно используется пользовательская норма, задаваемая общими параметрами интерфейса, а также, в случае с жестким решателем, ForwardDiff/DiffEqDiffTools для вычисления якобиана и базовые линейные разложения для нахождения линейного решения. Для вашего собственного типа вам, скорее всего, потребуется предоставить более оптимальную форму нормы, якобиана или вычислений для нахождения линейного решения, чтобы в полной мере использовать возможности параллелизма.

Пакеты GPUArrays.jl (CuArrays.jl), ArrayFire.jl и DistributedArrays.jl были протестированы и работают в различных формах, однако последний из них пока еще не рекомендуется для практического применения.

Следующий вопрос — имеют ли эти возможности значение. Как правило, эффект от параллелизма заметен лишь в больших системах. При использовании многопоточных массивов для трансляции мы ориентируемся на значение N>1000, хотя в руководстве Sundials указано значение N>100,000. Для методов Рунге-Кутта высокого порядка это значение, скорее всего, будет меньше, чем предлагаемое Sundials, так как в каждом внутреннем шаге больше операций. Однако, как всегда, для более точной оценки требуется проведение дополнительных тестов производительности, и значение будет зависеть от решаемой задачи. Чтобы применение GPU было целесообразным, как правило, требуется активное выполнение параллельных операций в пользовательской функции f, например перемножение матриц для вычисления шаблона в ДУЧП. Если вы просто решаете ОДУ поэлементно на большом массиве, скорее всего, особой пользы не будет либо скорость упадет из-за особенностей работы GPU. В случае с DistributedArrays применение параллельного режима для нахождения линейных решений должно быть действительно оправдано и поэтому рекомендуется только в случаях, когда данные задачи не помещаются в памяти или используется жесткий решатель для нахождения линейного решения методом Крылова.

Для перемножения и разложения матриц в Julia применяется базовая реализация BLAS. Эта библиотека автоматически работает в многопоточном режиме и ускоряет внутренние линейные алгебраические вычисления DifferentialEquations.jl. Однако для оптимального результата число используемых потоков BLAS должно соответствовать числу физических, а не логических ядер. Дополнительные сведения см. в описании этой проблемы.

Проверить число потоков BLAS можно так:

Чтобы напрямую задать 4 потока, используйте такое выражение:

Кроме того, иногда библиотека MKL от Intel может быть более быстрой реализацией BLAS, чем стандартная реализация BLAS, входящая в состав Julia (OpenBLAS). Чтобы переключиться на другую реализацию BLAS, можно воспользоваться пакетом MKL.jl, который ускоряет подпрограммы линейной алгебры. Это делается следующим образом:

Сначала проверьте наличие ошибок. Решатели проходят множество тестов на сходимость, поэтому если при их использовании возникает проблема, причина кроется либо в работе численных методов, либо в ошибке пользователя (последнее случается чаще, однако обратите внимание на раздел этой страницы, посвященный обычным числовым ошибкам). Самой частой причиной сообщений о медленном выполнении кода являются пользовательские ошибки в функции f, приводящие к расходимости решения.

Если у вас нет ошибок, отлично! Тогда можно использовать стандартные приемы оптимизации кода Julia. Некоторые советы и указания можно найти в руководстве по оптимизации кода DiffEq.

В первую очередь следует убедиться в том, что ваша функция не выделяет память. Если ваша система небольшая (<=100 ОДУ, СДУ, ДРУ или ДАУ), настройте для нее использование StaticArrays.jl. Это демонстрируется в руководстве по ОДУ со статическими матрицами. Статические векторы и массивы размещаются в стеке, поэтому новые массивы создаются без лишнего расходования памяти и компилятору не приходится размещать в куче временные данные (а это требует больше всего ресурсов). Для них предусмотрена специализированная, очень быстрая диспетчеризация для арифметических и других операций, таких как LU-разложения, поэтому по возможности лучше использовать их. Однако с чрезмерным увеличением размера их эффективность снижается.

Для больших данных следует использовать синтаксис на месте (in-place) f(du,u,p,t) и сделать так, чтобы функция не выделяла память. Если известно u0, можно сделать следующее.

В результате количество выделений памяти и объем выделяемой памяти должны быть близкими к нулю. Если это не так, возможно, имеет место нестабильность типов или есть временные массивы. Для выявления нестабильности типов следует выполнить следующий код:

После этого проверьте в выходных данных, есть ли типы, которые не выводятся компилятором, и исправьте их. Если есть глобальные переменные, их следует сделать константами (const). Что касается выделения памяти, обычно оно происходит по следующим причинам.

Пример оптимизации функции, получаемой в результате дискретизации ДУЧП, см. в этой записи блога.

Жесткие решатели требуют решения нелинейного уравнения на каждом шаге. Для этого требуется выполнить ряд шагов по методу Ньютона. По умолчанию эти методы предполагают, что якобиан является плотным, автоматически вычисляют якобиан и производят плотное разложение. Однако зачастую может быть желательно использовать альтернативные варианты, лучше подходящие для вашей задачи.

Во-первых, при наличии возможности рекомендуется передавать функцию для вычисления якобиана. Этот вопрос рассматривается в разделе о перегрузках производительности. Якобианы особенно полезны при использовании методов Розенброка.

Во-вторых, если ваш якобиан не плотный, не следует использовать плотный якобиан! Вместо этого, если применяется библиотека *DiffEq, следует указать линейный решатель и (или) jac_prototype для матричной формы, а для Sundials.jl следует изменить параметр linear_solver. Подробные сведения см. в разделе о решении ОДУ с помощью Sundials.

В настоящее время рекомендуемым методом для жестких задач с большими разреженными якобианами является QNDF. В качестве jac_prototype следует указать особую матрицу, например ленточную или тридиагональную, если она отвечает особой структуре. Если известно лишь то, что якобиан разреженный, с помощью автоматического определения разреженности можно установить шаблон разреженности. Дополнительные сведения см. в руководстве по жестким ОДУ. Наконец, LinSolveGMRES() может помочь в том случае, если для большой матрицы не удается получить шаблон разреженности или если из-за разреженности матрица не помещается в памяти. Как и ранее, справку по работе с дискретизациями ДУЧП можно найти в этой записи блога.

В этой записи на форуме подробно обсуждается, как быть с разрывами в функции ОДУ и как использовать дополнительную информацию для ускорения работы решателя.

Это делается несколькими способами. Самый простой из них — переключение между двумя функциями с помощью параметра. Например:

Затем в обратном вызове можно сделать так, чтобы функция affect! изменяла значение integrator.prob.p. Например, чтобы оно изменялось при u[2]<0.5:

Оно будет изменяться между двумя вариантами ОДУ для du1 в этот момент. Другой способ — привести функции ОДУ к одному типу посредством FunctionWrappers.jl, но в этом нет необходимости. В современных процессорах предусмотрено прогнозирование ветвления, поэтому условное выражение может свободно выполняться, если можно спрогнозировать, какая ветвь будет выбрана. В данном случае почти каждый вызов f идет по пути p==0 до обратного вызова, по достижении которого почти всегда выбирается путь else. Таким образом процессор эффективно избавляется от лишней вычислительной нагрузки, так что ваши усилия по оптимизации будут излишними (если только изменение не происходит на каждом шаге, но даже в таком случае это, скорее всего, будет наименее затратная часть вычислений).

Самые полезные параметры — это погрешности abstol и reltol. Они сообщают внутреннему механизму адаптивной дискретизации по времени, насколько точным должно быть решение. Как правило, reltol — это относительная точность, а abstol — точность при u, близком к нулю. Эти погрешности являются локальными и поэтому не гарантируются на глобальном уровне. Однако на практике общая точность решения обычно на 1—​2 знака ниже относительных погрешностей. Таким образом, для значений по умолчанию abstol=1e-6 и reltol=1e-3 можно ожидать общую точность на уровне примерно 1—​2 знаков. Это всеобщий стандарт, который применяется к собственным методам Julia, инкапсулированным методам Фортрана и C++, вызовам MATLAB, Python, R и т. д.

Да, причина в том, что численное решение ОДУ не является точным решением. Решить эту проблему можно несколькими путями. Один из них — получение более точного решения. Например, вместо

используйте

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

Другой путь — использование обратного вызова. В библиотеке обратных вызовов есть ряд готовых обратных вызовов, которые отвечают за такие задачи, как проецирование на многообразие и сохранение положительности.

Это невозможно. Для чисел с плавающей запятой не следует использовать погрешность менее abstol=1e-14 и reltol=1e-14. Если требуется меньшая погрешность, используйте числа произвольной точности, такие как BigFloats или ArbFloats.jl.

Да, они совместимы с автоматическим дифференцированием! Дополнительные сведения см. на странице, посвященной анализу чувствительности.

Если алгоритм не различает зависимые от параметров события, нужно просто использовать в качестве элементов начального условия дуальные числа. Если алгоритм использует дуальные числа, время также должно быть задано дуальными числами.

Чтобы продемонстрировать это на практике, предположим, что нужно найти якобиан решения для уравнения Лотки-Вольтерры при t=10 относительно параметров.

Эта функция принимает новые параметры и по завершении выдает решение. Начальное условие задано в виде eltype(p).([1.0,1.0]), чтобы применялись дуальные числа, если p — это массив чисел Dual. То же самое сделано для временного интервала, просто чтобы показать, что нужно делать, если есть зависимые от параметров события. После этого можно получить якобиан с помощью ForwardDiff.jl:

и сравнить его с FiniteDiff.jl:

Причина в том, что вы используете кэш, который несовместим с автоматическим дифференцированием посредством ForwardDiff.jl. Допустим, используется следующая функция ОДУ:

Здесь применяется временный кэшируемый массив, чтобы при перемножении матриц не выделялась память. При автоматическом дифференцировании типом элементов u являются числа Dual, поэтому A*u дает числа Dual и при попытке записи в tmp происходит ошибка. Есть два способа избежать этого. Первый способ проще и заключается в отключении автоматического дифференцирования с помощью параметра autodiff=false в решателе. Этот параметр есть у каждого решателя, использующего автоматическое дифференцирование. То есть решить задачу можно так:

Для вычисления якобианов будет использоваться резервный метод, а именно численное дифференцирование (DiffEqDiffTools.jl).

Для решения этой проблемы мы могли бы воспользоваться функциями get_tmp и dualcache из пакета PreallocationTools.jl, например:

Скорее всего, причина в том, что вы используете матрицу Якоби с меняющейся структурой разреженности, что несовместимо с линейным решателем по умолчанию для разреженных матриц. Если линейный решатель обнаружит эту проблему, появится сообщение об ошибке:

либо

Хотя ошибка Error: SingularException также возможна, если линейный решатель не обнаружил изменения структуры разреженности. Для решения этой проблемы необходимо отключить кэширование символического разложения, например:

Дополнительные сведения о возможных линейных решателях см. в документации LinearSolve.jl.

Отладчик несовместим с функцией llvmcall, используемой в форме AutoSpecialize для ускорения компиляции. Чтобы использовать отладчик, задействуйте форму FullSpecialize. То есть измените prob = ODEProblem(lorenz!,u0,tspan) на prob = ODEProblem{true, SciMLBase.FullSpecialize}(lorenz!,u0,tspan). Мы планируем исправить эту ошибку, но пока для любых случаев должно быть достаточно обходного решения.