//файл включения для векторного и матричного типов
#include "matrix.h"
/*для удобства переопределим шаблонные типы векторов и матриц для работы с числовыми объектами - вещественными числами двойной точности */
typedef vector<double> dvector;
typedef matrix<double> dmatrix;
/*Класс дифференциальных уравнений */
class DEqu
{
double t0, //начальное время
t1, //конечное время
step; //шаг по времени
dvector x0; //вектор начальных условий
/*Будем рассматривать класс уравнений, которые удается разрешить относительно старшей производной. В этом случае при представлении уравнения n-го порядка в виде системы уравнений первого порядка и решении этой системы нам нужна только одна функция - для правой части последнего уравнения в системе, вычисляющая его правую часть по уже вычисленным до этого значениям младших производных и самой функции. Для первого уравнения правая часть - в нулевой момент времени -это начальное значение первой производной, для второго уравнения - второй и т.д., а правая часть для последнего уравнения может включать зависимость от времени, значение функции и всех младших производных. Для унификации ввода этой функции в интерактивном режиме надо программировать интерпретатор. Так как нам этого делать не хотелось, то использовать настоящий класс сможет только программирующий пользователь - ему придется запрограммировать функцию правой части и указать ее имя при вызове соответствующего метода решения системы уравнений. Прототип этой функции объявлен через указатель на функцию, чтобы можно было использовать в качестве аргумента другой функции */
dvector (*dxbydt)(double t, dvector x);
public:
/*конструктор принимает в качестве параметров все вышеперечисленное */
/*Прототипы функций, реализующих методы Эйлера, Рунге-Кутты второго и четвертого порядка точности и Адамса. Функции возвращают матрицы с количеством столбцов, равным порядку системы, то есть в нулевом столбце - значения функции, а в последующих – значения производных соответствующего номеру столбца порядка */
dmatrix EilerSol();
dmatrix RK2Sol();
dmatrix RK4Sol();
dmatrix Adams4Sol();
/*набор функций для доступа и модификации параметров системы*/
/*Решение по методу Рунге-Кутта второго порядка (метод трапеций, или метод Хойна)*/
dmatrix DEqu::RK2Sol()
{
//определяем количество вычисляемых точек
long StepCount=(TEnd()-TStart())/StepByT();
//создаём матрицу для хранения решений
dmatrix x(StepCount+1, 1);
//заносим начальные условия в 0-й вектор матрицы
x[0]=XStart();
for(long i=0;i<StepCount;i++)
{
/*вычисляем вектор-коэффициенты следующей точки решения*/
dvector k1=dxbydt(TStart()+i*StepByT(),x[i]);
dvector k2=dxbydt(TStart()+(i+1)*StepByT(),
x[i]+StepByT()*k1);
//вычисляем вектор решения на следующем шаге
x[i+1]=x[i]+(StepByT()/2)*(k1+k2);
}
return x;//возвращаем матрицу-результат
}
/*Решение по методу Рунге-Кутта 4-го порядка точности*/
dmatrix DEqu::RK4Sol()
{
//количество точек решения
long StepCount=(TEnd()-TStart())/StepByT();
dmatrix x(StepCount+1, 1);//матрица решения
x[0]=XStart();//начальные условия
for(long i=0;i<StepCount;i++)
{
//вычисляем вектор-коэффициенты решения
dvector k1=dxbydt(TStart()+i*StepByT(),x[i]);
dvector k2=dxbydt(TStart()+(i+0.5)*StepByT(),
x[i]+(StepByT()/2)*k1);
dvector k3=dxbydt(TStart()+(i+0.5)*StepByT(),
x[i]+(StepByT()/4)*(k1+k2));
dvector k4=dxbydt(TStart()+(i+1)*StepByT(),
x[i]+StepByT()*(-k2+2*k3));
x[i+1]=x[i]+(StepByT()/6)*(k1+4*k3+k4); //следующая точка
}
return x;//возвращаем матрицу решений
}
//Метод Адамса
dmatrix DEqu::Adams4Sol()
{
//определяем количество точек решения
long StepCount=(TEnd()-TStart())/StepByT();
dmatrix x(StepCount+1,1);/*матрица решения - массив векторов*/
/*преобразуем указатель на данный объект в указатель на объект для решения системы дифференциальных уравнений методом Рунге-Кутта четвёртого порядка и получаем объект соответствующего типа вызовом конструктора копирования */
DEqu obj=*this;
/*если количество точек решения таково, что данную систему нельзя решить методом Адамса-Башфорта, возвращаем решение по методу Рунге-Кутта*/
if(StepCount<=3) return obj.RK2Sol();
/*В противном случае модифицируем параметры системы (конечное время) так, чтобы по методу Рунге-Кутта четвёртого порядка находились только первые четыре решения */
obj.TEnd()=TStart()+3*StepByT();
dmatrix tempx=obj.RK4Sol();/*получаем первые четыре точки*/
for(long i=0;i<4;i++)/*и переписываем их в векторы результирующей матрицы*/
x[i]=tempx[i];
/*начиная с пятой точки, работает собственно метод Адамса*/
for(long i=3;i<StepCount;i++)
x[i+1]=x[i]+(StepByT()/24)*
(
55*dxbydt(TStart()+i*StepByT(),x[i])-
59*dxbydt(TStart()+(i-1)*StepByT(),x[i-1])+
37*dxbydt(TStart()+(i-2)*StepByT(),x[i-2])-
9*dxbydt(TStart()+(i-3)*StepByT(),x[i-3])
);
return x;//возвращаем результирующую матрицу
}
//Тестовые функции
dvector testfunc(double x, dvector y)
{
//Решение уравнения Бесселя
dvector dy=2;
dy[0]=y[1];/*начальное условие для первой производной*/
