1

C/c++ をやってからしばらく経ちましたが、次の ODE セットを解くために gsl ライブラリを使用して 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 表記では my y[0]=u, y[1]==up となり、上記の RHS は f[0] と f[1] を定義します。これらの定義から、ヤコビアンと dfdr を計算できます (通常、それらの「時間」変数は「r」ではなく「t」と呼ばれます)。これを行う理由は、Mathematica の速度に問題があるためです。私は、ODE ソルバーに関するドキュメントの最後にある gsl サンプル コードを取り上げ、次のように私の問題に適用しようとしました。

 #include <stdio.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_matrix.h>
 #include <gsl/gsl_odeiv2.h>
 #include <complex.h>


 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;
 }

これは数値を与えていますが、WorkingPrecision と PrecisionGoal が低い場合でも、Mathematica NDSolve が与える数値と同じではありません。私が行ったことに誤りはありますか？

4

1 に答える 1

4

括弧が足りないだけだと思い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ます。

于 2012-10-24T17:21:46.773 に答える