Обзор доступных библиотек для численного решения жёстких ОДУ

Создавая дополнения к отечественной математической программе SMath Studio ^[1], я нашёл в сети ряд библиотек, которые можно было бы использовать в своих программах. Предлагаю небольшой их обзор.

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

Большинство представленных ниже функций, если не оговорено особо, можно привести к одному формату вызова (по аналогии с Mathcad):

ode_solver( init, x1, x2, intvls, D(t, x) )

где:

• init — вектор начальных условий,
• (x1, x2) — отрезок интегрирования,
• intvls — количество интервалов на отрезке,
• D(t, x) — система ОДУ.

1. Intel ODE Solvers Library

Содержит следующие функции: rkm9st(), mk52lfn(), mk52lfa(), rkm9mkn(), rkm9mka().

• rkm9st() — a specialized routine for solving non-stiff and middle-stiff ODE systems using the explicit method, which is based on the 4th order Merson’s method and the 1st order multistage method of up to and including 9 stages with stability control.
• mk52lfn() — a specialized routine for solving stiff ODE systems using the implicit method based on L-stable (5,2)-method with the numerical Jacobi matrix, which is computed by the routine.
• mk52lfa() — a specialized routine for solving stiff ODE systems using the implicit method based on L-stable (5,2)-method with numerical or analytical computation of the Jacobi matrix. The user must provide a routine for this computation.
• rkm9mkn() — a specialized routine for solving ODE systems with a variable or a priori unknown stiffness; automatically chooses the explicit or implicit scheme in every step and computes the numerical Jacobi matrix when necessary.
• rkm9mka() — a specialized routine for solving ODE systems with a variable or a priori unknown stiffness; automatically chooses the explicit or implicit scheme in every step. The user must provide a routine for numerical or analytical computation of the Jacobi matrix.

Библиотека написана на C со всеми вытекающими отсюда зависимостями. Доступны 32- и 64-разрядные версии библиотеки (libiode_ia32.lib и libiode_intel64.lib).

/******************************************************************************* 
!                              INTEL CONFIDENTIAL 
!   Copyright(C) 2007-2008 Intel Corporation. All Rights Reserved. 
!   The source code contained  or  described herein and all documents related to 
!   the source code ("Material") are owned by Intel Corporation or its suppliers 
!   or licensors.  Title to the  Material remains with  Intel Corporation or its 
!   suppliers and licensors. The Material contains trade secrets and proprietary 
!   and  confidential  information of  Intel or its suppliers and licensors. The 
!   Material  is  protected  by  worldwide  copyright  and trade secret laws and 
!   treaty  provisions. No part of the Material may be used, copied, reproduced, 
!   modified, published, uploaded, posted, transmitted, distributed or disclosed 
!   in any way without Intel's prior express written permission. 
!   No license  under any  patent, copyright, trade secret or other intellectual 
!   property right is granted to or conferred upon you by disclosure or delivery 
!   of the Materials,  either expressly, by implication, inducement, estoppel or 
!   otherwise.  Any  license  under  such  intellectual property  rights must be 
!   express and approved by Intel in writing.
! 
!******************************************************************************
!   
!  Header file for Intel(R) ODE Solvers
! 
!*******************************************************************************/

#ifndef _INTEL_ODE_H_
#define _INTEL_ODE_H_

#ifdef __cplusplus
extern "C" {
#endif /* __cplusplus */

void dodesol(int*,int*,double*,double*,double*,void*,void*,
			 double*,double*,double*,double*,double*,int*,int*);
void dodesol_rkm9st(int*,int*,double*,double*,double*,void*,
					double*,double*,double*,double*,double*,int*);
void dodesol_mk52lfn(int*,int*,double*,double*,double*,void*,
					 double*,double*,double*,double*,double*,int*,int*);
void dodesol_mk52lfa(int*,int*,double*,double*,double*,void*,void*,
					 double*,double*,double*,double*,double*,int*,int*);
void dodesol_rkm9mkn(int*,int*,double*,double*,double*,void*,
					 double*,double*,double*,double*,double*,int*,int*);
void dodesol_rkm9mka(int*,int*,double*,double*,double*,void*,void*,
					 double*,double*,double*,double*,double*,int*,int*);

#ifdef __cplusplus
}
#endif /* __cplusplus */

#endif /* _INTEL_ODE_H_ */

Дополнение ODE Solvers ^[2] демонстрирует работу с этой библиотекой из c# кода.

Ссылки:
1. Intel® Ordinary Differential Equations Solver Library ^[3].
2. Исходники дополнения ODESolvers ^[4].

2. GNU Scientific Library (GSL)

Содержит следующие функции: gslrk2(), gslrk4(), gslrkf45(), gslrkck(), gslrk8pd(), rk1imp(), rk2imp(), rk4imp(), bsimp(), msadams(), msbdf().

Часть из них требует дополнительные параметры для работы (Якобиан). Те, которые мне удалось привести к общему виду:

Solvers for Non-Stiff Systems:

• gslrk2() — explicit embedded Runge-Kutta (2, 3) method.
• gslrk4() — explicit 4th order (classical) Runge-Kutta. Error estimation is carried out by the step doubling method.
• gslrkf45() — explicit embedded Runge-Kutta-Fehlberg (4, 5) method.
• gslrkck() — explicit embedded Runge-Kutta Cash-Karp (4, 5) method.
• gslrk8pd() — explicit embedded Runge-Kutta Prince-Dormand (8, 9) method.

Для работы с функциями используется универсальный интерфейс ^[5], где конкретный тип решателя задаёт шаговую функцию. Дополнение GNUScientificLibrary ^[6] демонстрирует работу с этой библиотекой из c# кода.

Не так просто сделать сборку библиотеки под Windows. Я использовал инструкцию с одного сайта, который сейчас недоступен. Тем не менее, в репозитории дополнения вы сможете найти 32- и 64-разрядные версии dll для GSL 1.16. Там находится вся библиотека, а не только решатели ОДУ.

Ссылки:
1. GSL. Ordinary Differential Equations ^[7].
2. Исходники дополнения GNUScientificLibrary ^[8].

3. Matlab C++ Math Library 2.1 (Win32)

Да, вы можете использовать эту старую версию run-time библиотек для расчётов. Более того, она может быть установлена по относительным путям, т.е. можно просто положить содержимое оригинального дистрибутива (~28 Мб в развёрнутом виде) рядом со своей программой. Правда при вызове функций придётся использовать SetCurrentDirectory() с прямым указанием на место расположения «binwin32». Я так делаю в своём дополнении.

Содержит следующие функции: ode23(), ode45(), ode113(), ode15s(), ode23s().

• ode23() — solve nonstiff differential equations; low order method,
• ode45() — solve nonstiff differential equations; medium order method,
• ode113() — solve nonstiff differential equations; variable order method,
• ode15s() — solve stiff differential equations and DAEs; variable order method,
• ode23s() — solve stiff differential equations; low order method.

Дополнение MatlabCppMathLibrary ^[9] демонстрирует работу с этой библиотекой из c# кода.

Ссылки:
1. Ordinary Differential Equations ^[10].
2. MATLAB C++ Math Library. User's Guide. Version 2.1 ^[11] (pdf).
3. MATLAB C++ Math Library. Reference. Version 2 ^[12] (pdf).
4. Исходники дополнения MatlabCppMathLibrary ^[13].

4. Octave C++ Math Library (Win32)

Примерно то же самое, что и Matlab C++ Math Library, но со своими тараканами. К сожалению, работу с этой библиотекой я одолел только частично. Дополнение OctaveCppMathLibrary ^[14] демонстрирует работу с этой библиотекой из c# кода.

Ссылки:
1. Ordinary Differential Equations ^[15].
2. Исходники дополнения OctaveCppMathLibrary ^[16].

5. DotNumerics

Содержит следующие функции: AdamsMoulton(), ExplicitRK45(), ImplicitRK5(), GearsBDF(). Эта библиотека портирована для .Net с фортрана. Она понравилась мне больше всего. Работает достаточно быстро.

Solvers for Non-Stiff Systems:

• AdamsMoulton() — solves an initial-value problem for nonstiff ordinary differential equations using the Adams-Moulton method.
• ExplicitRK45() — solves an initial-value problem for nonstiff ordinary differential equations using the explicit Runge-Kutta method of order (4)5.

Solvers for Stiff Systems:

• ImplicitRK5() — solves an initial-value problem for stiff ordinary differential equations using the implicit Runge-Kutta method of order 5.
• GearsBDF() — solves an initial-value problem for stiff ordinary differential equations using the Gear’s BDF method.

Имеется много перегрузок для различных форматов вызова функций. Дополнение DotNumerics ^[17] демонстрирует работу с этой библиотекой.

Ссылки:
1. DotNumerics ^[18].
2. Исходники дополнения DotNumerics ^[19].

Как выглядит модель амплитудного детектора в SMath Studio при решении ОДУ с помощью функции GearsBDF():

SMath Studio

Автор: ViacheslavMezentsev

Источник ^[20]