У меня есть ODE типа:
x'(t) = a(t)x+g(t)
Который я пытаюсь решить. Единственный Пример GSL ODE не очень полезно, потому что единственный коэффициент (\ mu) не зависит от времени.
Этот вопрос был ответил в списке рассылки GSL однако ответ очень неясен — g(t)
игнорируется, и не было объяснено, как включить a(t)
в func
(должно быть передано в *params
?).
Есть ли какой-нибудь пример, который я могу увидеть, где такое ODE решается с помощью GSL?
ОБНОВИТЬ: Как указывалось ниже, на это был дан ответ в списке рассылки GSL. Вот полный пример программы, как это делается:
#include <stdio.h>
#include <math.h>
#include "gsl/gsl_errno.h"#include "gsl/gsl_matrix.h"#include "gsl/gsl_odeiv2.h"int func(double t, const double y[], double f[], void *params) {
f[0] = -t* y[0];
return GSL_SUCCESS;
}
int jac(double t, const double y[], double *dfdy, double dfdt[], void
*params) {
gsl_matrix_view dfdy_mat = gsl_matrix_view_array(dfdy, 1, 1);
gsl_matrix * m = &dfdy_mat.matrix;
gsl_matrix_set(m, 0, 0, -t);
dfdt[0] = -1;
return GSL_SUCCESS;
}
int main(void) {
double mu = 0;
gsl_odeiv2_system sys = { func, jac, 1, &mu };
gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new(&sys,
gsl_odeiv2_step_rk1imp, 1e-7, 1e-7, 0.0);
int i;
double t = 0.0, t1 = 2.0;
const double x0 = 1.0;
double y[1] = {x0};
const int N = 100;
printf("time\t \tapprox solution \t exact solution\n");
for (i = 0; i <= N; i++) {
double ti = i * (t1 / N);
int status = gsl_odeiv2_driver_apply(d, &t, ti, y);
if (status != GSL_SUCCESS) {
printf("error, return value=%d\n", status);
break;
}
printf("%.5e\t%.5e\t\t%.5e\n", t, y[0], x0*exp(-0.5*t*t));
}
gsl_odeiv2_driver_free(d);
printf("\n...and done!\n");
return 0;
}
На этот вопрос уже ответили здесь http://lists.gnu.org/archive/html/help-gsl/2012-06/msg00036.html
Если вы не ограничены GSL и / или C, вы можете использовать http://odeint.com — современная библиотека C ++ для ODE. Odeint является частью поддержки, поэтому он может быть уже установлен в вашей системе или может быть легко установлен в большинстве менеджеров пакетов для дистрибутивов Linux.
Вы можете просто определить свои коэффициенты и ODE и использовать, например, метод RK4:
double coef_a( double t ) { /* return a(t) */ };
double coef_g( double t ) { /* return b(t) */ };
typedef std::array< double , 1 > state_type;
double ode( state_type const &x , state_type &dxdt , double )
{
dxdt[0] = coef_a( t ) * x[0] + coef_g( t );
}
state_type x;
double t_state , t_end , dt;
// initialize x
integrate_const( runge_kutta< state_type >() , ode , x , t_start , dt , t_end );