Создавая дополнения к отечественной математической программе SMath Studio, я нашёл в сети ряд библиотек, которые можно было бы использовать в своих программах. Предлагаю небольшой их обзор.
Стандартный 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 демонстрирует работу с этой библиотекой из c# кода.
Ссылки:
1. Intel® Ordinary Differential Equations Solver Library.
2. Исходники дополнения ODESolvers.
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.
Для работы с функциями используется универсальный интерфейс, где конкретный тип решателя задаёт шаговую функцию. Дополнение GNUScientificLibrary демонстрирует работу с этой библиотекой из c# кода.
Не так просто сделать сборку библиотеки под Windows. Я использовал инструкцию с одного сайта, который сейчас недоступен. Тем не менее, в репозитории дополнения вы сможете найти 32- и 64-разрядные версии dll для GSL 1.16. Там находится вся библиотека, а не только решатели ОДУ.
Ссылки:
1. GSL. Ordinary Differential Equations.
2. Исходники дополнения GNUScientificLibrary.
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 демонстрирует работу с этой библиотекой из c# кода.
Ссылки:
1. Ordinary Differential Equations.
2. MATLAB C++ Math Library. User's Guide. Version 2.1 (pdf).
3. MATLAB C++ Math Library. Reference. Version 2 (pdf).
4. Исходники дополнения MatlabCppMathLibrary.
4. Octave C++ Math Library (Win32)
Примерно то же самое, что и Matlab C++ Math Library, но со своими тараканами. К сожалению, работу с этой библиотекой я одолел только частично. Дополнение OctaveCppMathLibrary демонстрирует работу с этой библиотекой из c# кода.
Ссылки:
1. Ordinary Differential Equations.
2. Исходники дополнения OctaveCppMathLibrary.
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 демонстрирует работу с этой библиотекой.
Ссылки:
1. DotNumerics.
2. Исходники дополнения DotNumerics.
Как выглядит модель амплитудного детектора в SMath Studio при решении ОДУ с помощью функции GearsBDF():
Автор: ViacheslavMezentsev
