odeint неявный простой пример Эйлера
Я заинтересован в решении системы ODE с библиотекой odeint, используя неявную схему, и у меня есть трудности с реализацией простого implicit_euler пример.
Глядя на документацию, мне удалось сделать работу явным степпером, адаптивным, а также rosenbrock4 шаговый. Первые кажутся полуявными. Поэтому я был заинтересован в реализации полностью неявной схемы (и в то же время извлекать матрицу Якоби на каждом временном шаге). Но мне не удалось найти документацию и рабочие примеры этого степпера. Что у меня есть
typedef boost::numeric::ublas::vector< double > vector_type;
typedef boost::numeric::ublas::matrix< double > matrix_type;`
struct stiff_system
{
void operator()( const vector_type &x , vector_type &dxdt , double /* t */ )
{
dxdt[ 0 ] = -101.0 * x[ 0 ] - 100.0 * x[ 1 ];
dxdt[ 1 ] = x[ 0 ];
}
};
struct stiff_system_jacobi
{
void operator()( const vector_type & /* x */ , matrix_type &J , const double & /* t */ , vector_type &dfdt )
{
J( 0 , 0 ) = -101.0;
J( 0 , 1 ) = -100.0;
J( 1 , 0 ) = 1.0;
J( 1 , 1 ) = 0.0;
dfdt[0] = 0.0;
dfdt[1] = 0.0;
}
};
typedef implicit_euler< double > stepper_IE;
vector_type inout( 2 , 1.0 );
size_t steps = integrate_const( stepper_IE() ,
std::make_pair( stiff_system() , stiff_system_jacobi() ) ,
inout , 0.0 , 5.0 , 0.01, streaming_observer( std::cout, x_vec , times ));
Ошибка заключается в следующем:
C:\boost_1_55_0\boost\numeric\odeint\stepper\implicit_euler.hpp:94: ошибка: C2064: термин не оценивает функцию, принимающую 3 аргумента
класс не определяетoperator()‘или определенный пользователем оператор преобразования в указатель на функцию или ссылку на функцию, который принимает соответствующее количество аргументов
Мой вопрос сейчас: кто-то знает, как заставить это работать, или может кто-то указать мне на более подробную документацию, чем эта:
или этот:
Спасибо
Решение
К сожалению, неявный метод Эйлера и решатель Розенброка (еще один неявный решатель) не имеют одинакового интерфейса. Подробно, неявный Эйлер ожидает функцию для якобиана с этой сигнатурой
void jacobian( const state_type &x , matrix_type &jacobi , const value_type t );
Следовательно, вам нужно изменить свое определение stiff_system_jacobi в
struct stiff_system_jacobi
{
void operator()( const vector_type & , matrix_type &J , const double & ) const
{
J( 0 , 0 ) = -101.0;
J( 0 , 1 ) = -100.0;
J( 1 , 0 ) = 1.0;
J( 1 , 1 ) = 0.0;
}
};
Если ваша система действительно не автономна, вам нужно улучшить свой тип состояния с помощью одной дополнительной координаты, которая представляет время и имеет тривиальную динамику dt / dt = 1.
Другие решения