使用gsl编写Runge-Kutta ODE求解器

自从我做了任何C / c ++以来已经有一段时间了，但我想用gsl库编写一个ODE求解器来解决下面的ODE集

$$ u'(r)=up(r)$$ $$ up'(r)=-(2*(r-1)/(r*(r-2)))*up(r)-((r*r/((r-2)*(r-2)))-(2/r*(r-2)))*u(r) $$ 

所以在gsl表示法中我的y [0] = u，y [1] == up，并且上面的RHS定义了f [0]和f [1]。 然后可以从这些定义中计算雅可比行列式和dfdr（通常它们的“时间”变量称为“t”而不是“r”）。 这样做的原因是因为我对Mathematica有速度问题。 我在ODE解算器的文档末尾使用了gsl示例代码，并尝试将其调整到我的问题，如下所示：

  #include  #include  #include  #include  #include  int func (double r, const double y[], double f[], void *params) { double mu = *(double *)params; f[0] = y[1]; f[1] = -(2*(r-1)/(r*(r-2)))*y[1]-((r*r/((r-2)*(r-2)))-(2/r*(r-2)))*y[0]; return GSL_SUCCESS; } /* void tester (double r) { double outer=-((r*r/((r-2)*(r-2)))-(2/(r*(r-2)))); printf ("%.5e \n", outer); } */ int jac (double r, const double y[], double *dfdy, double dfdt[], void *params) { double mu = *(double *)params; gsl_matrix_view dfdy_mat = gsl_matrix_view_array (dfdy, 2, 2); gsl_matrix * m = &dfdy_mat.matrix; gsl_matrix_set (m, 0, 0, 0.0); gsl_matrix_set (m, 0, 1, 1.0); gsl_matrix_set (m, 1, 0,-((r*r/((r-2)*(r-2)))-(2/(r*(r-2))))); gsl_matrix_set (m, 1, 1, -(2*(r-1)/(r*(r-2)))); dfdt[0] = 0.0; dfdt[1] =((1/(r*r))+(1/((r-2)*(r-2))))*y[1]-((4*(1-r)/(r*r*(r-2)*(r-2)))+(4*r/((r-2)*(r-2)*(r-2))))*y[0]; return GSL_SUCCESS; } int main (void) { /* tester(5);*/ double om = 2; gsl_odeiv2_system sys = {func, jac, 2, &om}; gsl_odeiv2_driver * d = gsl_odeiv2_driver_alloc_y_new (&sys, gsl_odeiv2_step_rk8pd, 1e-6, 1e-6, 0.0); int i; double r = 10, r1 = 100; double y[2] = { 0.0000341936, -0.0000572397 }; for (i = 1; i <= 90; i++) { double ri = 10 + i; int status = gsl_odeiv2_driver_apply (d, &r, ri, y); if (status != GSL_SUCCESS) { printf ("error, return value=%d\n", status); break; } printf ("%.5e %.5e %.5e\n", r, y[0], y[1]); } gsl_odeiv2_driver_free (d); return 0; } 

这是给出数字，但它们与Mathematica NDSolve给出的数字不同，即使我的WorkingPrecision和PrecisionGoal较低。 我做了什么错误？

Answer 1

你只是在r*(r-2)周围缺少括号

 ((r*r/((r-2)*(r-2)))-(2/r*(r-2)))*y[0]; 

（也在up'(r)=-(2*(r-1)/(r*(r-2)))*up(r)-((r*r/((r-2)*(r-2)))-(2/r*(r-2)))*u(r) ），你在jac有相应的括号。