Программа предназначена для решения систем ОДУ, содержащих некоторую неопределённость или неоднозначность в начальных условиях или параметрах (т.е. в тех случаях, когда начальные условия задаются интервалами). Реализация метода основана на символьных вычислениях и использовании модели Тейлора.
Работа базируется главным образом на статье Markus Neher «From Interval Analysis to Taylor Models - An Overview». Математические операции над рядами описаны в статье N. Revol, K. Makino, M. Berz «Taylor Models and Floating-Point Arithmetic: Proof that Arithmetic Operations are Validated in COSY». Особое внимание уделено алгоритму перемножения рядов, которое описано в статье Martin Berz «Algorithms for Higher Order Automatic Differentiation in Many Variables with Applications to Beam Physics».
Программа является ресурсозатратной, так как она, во-первых, оперирует рядами и, во-вторых, дополнительно вычисляет вспомогательные таблицы коэффициентов для быстрого перемножения рядов. Потому программа писалась с изначальным расчётом на то, что она будет решать одну-единственную систему ОДУ за раз (смиритесь =3 ).
Разберём, как записать систему, описанную в статье «From Interval Analysis to Taylor Models - An Overview».
Прежде всего следует понять, сколько здесь переменных, а сколько параметров. Подставив a и b в начальные условия, можно увидеть, что у системы нет параметров, есть только переменные. Записываем исходные интервалы в вектор (порядок важен!).
vector<interval<double> > initPoint;
initPoint.push_back(interval<double>(0.95, 1.05));
initPoint.push_back(interval<double>(-1.05, -0.95));Для инициализации системы требуется указать количество переменных и параметров системы и порядок ряда, в котором будут производиться расчёты. А также начальные условия системы, которые предаются в виде вектора интервалов.
equation<double> odu(2, 0, 18);
odu.initialFlow(&initPoint);Затем необходимо определить сами рассчитываемые функции. В классе equation объявляем все функции системы и помещаем их в вектор функций.
powerSeries<T> pFun1(vector<powerSeries<T> >&);
powerSeries<T> pFun2(vector<powerSeries<T> >&);
vector<mfunction> pFun = { &equation<T>::pFun1, &equation<T>::pFun2 };Определяем функции.
template <typename T>
powerSeries<T> equation<T>::pFun1(vector<powerSeries<T> > &v) {
return v[1];
}
template <typename T>
powerSeries<T> equation<T>::pFun2(vector<powerSeries<T> > &v) {
return v[0] * v[0];
}Наконец, осуществляем расчёт методом Рунге-Кутты 4-го порядка.
odu.RungeKutta(0, 6, 0.01);
Все последующие примеры взяты из статьи А.Ю. Морозова и Д.Л. Ревизникова «Алгоритм адаптивной интерполяции на основе kd-дерева для численного интегрирования систем ОДУ с интервальными начальными условиями»
Записываем начальные условия, сначала идут переменные и только затем параметры. Это важно!
vector<interval<double> > initPoint;
initPoint.push_back(interval<double>(0.9, 1.1));
initPoint.push_back(interval<double>(1.9, 2.1));
initPoint.push_back(interval<double>(0.7, 0.75));
equation<double> odu(2, 1, 4);
odu.initialFlow(&initPoint);Затем определяем функции.
Обратите внимание, что к переменным мы обращаемся через аргумент функции v[i], а вот к параметру системы через u[j]. Как определить j? У него тот же индекс, что и в векторе initPoint.
template <typename T>
powerSeries<T> equation<T>::pFun1(vector<powerSeries<T> > &v) {
return v[0] * (-0.9) + v[0] * v[1] * 0.5;
}
template <typename T>
powerSeries<T> equation<T>::pFun2(vector<powerSeries<T> > &v) {
return u[2] * v[1] + v[0] * v[1] * (-0.8);
}Рассчитываем систему с выводом графика с шагом печати 30*h.
odu.RungeKutta(0, 100, 0.01, true, 30);
На графике отражена проекция множества решений из трёхмерного пространства в двухмерное.
Немного о тригонометрии. Она пока не реализована, придётся раскладывать функции в ряд самостоятельно.
Определяем функции.
template <typename T>
powerSeries<T> equation<T>::pFun1(vector<powerSeries<T> > &v) {
return v[1];
}
template <typename T>
powerSeries<T> equation<T>::pFun2(vector<powerSeries<T> > &v) {
powerSeries<T> p3 = v[0] * v[0] * v[0];
powerSeries<T> p5 = p3 *v[0] * v[0];
powerSeries<T> p7 = p5*v[0] * v[0];
powerSeries<T> p9 = p7*v[0] * v[0];
return (v[0] - p3 / 6 + p5 / 120 + p7 / 5040)*(-1);
}И производим расчёт.
vector<interval<double> > initPoint;
initPoint.push_back(interval<double>(-1.0, 1.0));
initPoint.push_back(interval<double>(0, 1.0));
equation<double> odu(2, 0, 8);
odu.initialFlow(&initPoint);
odu.RungeKutta(0, 15, 0.01, true, 30);
equation(int nvar, int param, int order) – инициализация класса.
nvar – количество переменных системы;
param – количество параметров системы;
order – порядок ряда, в котором будут проводиться все последующие вычисления.
void initialFlow(vector<interval > *points) – инициализация начальных условий системы. В points должны находиться начальные интервалы переменных, затем параметров.
void RungeKutta(double tStart, double tEnd, double h, bool plot = false, int plotStep = 0, std::string filename = "function.dat") – расчёт системы методом Рунге-Кутты 4-го порядка.
tStart – начальное время расчёта системы;
tEnd – конечное время расчёта системы;
h – шаг по времени;
plot – указывается true, если при расчёте необходимо вывести значения системы в файл для дальнейшего построения графика. По умолчанию равно false;
plotStep – шаг печати, по умолчанию равен (1.0 / 2h);
filename – имя файла, в который будет выводиться значения системы во время расчётов. По умолчанию "function.dat".
void printPlot(std::string filename) – выведет в файл с именем filename состояние системы на текущий момент.
В отличие от класса equation, экземпляров класса multSerCoef можно создать неограниченное количество.
int order() – возвращает порядок ряда.
int realVariable() – возвращает количество переменных системы.
int variableEven() – возвращает количество переменных, нужных для построения таблицы коэффициентов. По факту это сумма числа переменных и параметров системы ОДУ, доведённых до чётного числа.
int realParameter() - возвращает количество параметров системы.
int serieSize() – возвращает количество членов, составляющих ряд.
void printTableC() и void printTableD() – выведет в консоль таблицы коэффициентов (см. соответствующую статью)