test.c

该文件为c++的数学函数库!是一个非常有用的编程工具.它含有各种数学函数,为科学计算、工程应用等程序编写提供方便！

/* ode-initval/test_odeiv.c
 * 
 * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
 * 
 * This program is free software; you can redistribute it and/or modify
 * it under the terms of the GNU General Public License as published by
 * the Free Software Foundation; either version 2 of the License, or (at
 * your option) any later version.
 * 
 * This program is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
 * General Public License for more details.
 * 
 * You should have received a copy of the GNU General Public License
 * along with this program; if not, write to the Free Software
 * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
 */

/* Author:  G. Jungman
 */
#include <config.h>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_test.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_ieee_utils.h>
#include <gsl/gsl_odeiv.h>

int rhs_linear (double t, const double y[], double f[], void *params);
int jac_linear (double t, const double y[], double *dfdy, double dfdt[],
                void *params);
int rhs_sin (double t, const double y[], double f[], void *params);
int jac_sin (double t, const double y[], double *dfdy, double dfdt[],
             void *params);
int rhs_exp (double t, const double y[], double f[], void *params);
int jac_exp (double t, const double y[], double *dfdy, double dfdt[],
             void *params);
int rhs_stiff (double t, const double y[], double f[], void *params);
int jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
               void *params);
void test_stepper_linear (const gsl_odeiv_step_type * T, double h,
                         double base_prec);
void test_stepper_sin (const gsl_odeiv_step_type * T, double h, double base_prec);
void test_stepper_exp (const gsl_odeiv_step_type * T, double h, double base_prec);
void test_stepper_stiff (const gsl_odeiv_step_type * T, double h, double base_prec);
void test_evolve_system_flat (gsl_odeiv_step * step,
                              const gsl_odeiv_system * sys,
                              double t0, double t1, double hstart,
                              double y[], double yfin[],
                              double err_target, const char *desc);

void test_evolve_system (const gsl_odeiv_step_type * T,
                         const gsl_odeiv_system * sys,
                         double t0, double t1, double hstart,
                         double y[], double yfin[],
                         double err_target, const char *desc);

void test_evolve_exp (const gsl_odeiv_step_type * T, double h, double err);
void test_evolve_sin (const gsl_odeiv_step_type * T, double h, double err);
void test_evolve_stiff1 (const gsl_odeiv_step_type * T, double h, double err);
void test_evolve_stiff5 (const gsl_odeiv_step_type * T, double h, double err);

/* RHS for a + b t */

int
rhs_linear (double t, const double y[], double f[], void *params)
{
  f[0] = 0.0;
  f[1] = y[0];
  return GSL_SUCCESS;
}

int
jac_linear (double t, const double y[], double *dfdy, double dfdt[],
            void *params)
{
  gsl_matrix dfdy_mat;
  dfdy_mat.size1 = 2;
  dfdy_mat.size2 = 2;
  dfdy_mat.tda = 2;
  dfdy_mat.data = dfdy;
  dfdy_mat.block = 0;
  gsl_matrix_set (&dfdy_mat, 0, 0, 0.0);
  gsl_matrix_set (&dfdy_mat, 0, 1, 0.0);
  gsl_matrix_set (&dfdy_mat, 1, 0, 1.0);
  gsl_matrix_set (&dfdy_mat, 1, 1, 0.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_lin = {
  rhs_linear,
  jac_linear,
  2,
  0
};


/* RHS for sin(t),cos(t) */

int
rhs_sin (double t, const double y[], double f[], void *params)
{
  f[0] = -y[1];
  f[1] = y[0];
  return GSL_SUCCESS;
}

int
jac_sin (double t, const double y[], double *dfdy, double dfdt[],
         void *params)
{
  gsl_matrix dfdy_mat;
  dfdy_mat.data = dfdy;
  dfdy_mat.size1 = 2;
  dfdy_mat.size2 = 2;
  dfdy_mat.tda = 2;
  dfdy_mat.block = 0;
  gsl_matrix_set (&dfdy_mat, 0, 0, 0.0);
  gsl_matrix_set (&dfdy_mat, 0, 1, -1.0);
  gsl_matrix_set (&dfdy_mat, 1, 0, 1.0);
  gsl_matrix_set (&dfdy_mat, 1, 1, 0.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_sin = {
  rhs_sin,
  jac_sin,
  2,
  0
};


/* RHS for a exp(t)+ b exp(-t) */

int
rhs_exp (double t, const double y[], double f[], void *params)
{
  f[0] = y[1];
  f[1] = y[0];
  return GSL_SUCCESS;
}

int
jac_exp (double t, const double y[], double *dfdy, double dfdt[],
         void *params)
{
  gsl_matrix dfdy_mat;
  dfdy_mat.data = dfdy;
  dfdy_mat.size1 = 2;
  dfdy_mat.size2 = 2;
  dfdy_mat.tda = 2;
  dfdy_mat.block = 0;
  gsl_matrix_set (&dfdy_mat, 0, 0, 0.0);
  gsl_matrix_set (&dfdy_mat, 0, 1, 1.0);
  gsl_matrix_set (&dfdy_mat, 1, 0, 1.0);
  gsl_matrix_set (&dfdy_mat, 1, 1, 0.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_exp = {
  rhs_exp,
  jac_exp,
  2,
  0
};


/* RHS for stiff example */

int
rhs_stiff (double t, const double y[], double f[], void *params)
{
  f[0] = 998.0 * y[0] + 1998.0 * y[1];
  f[1] = -999.0 * y[0] - 1999.0 * y[1];
  return GSL_SUCCESS;
}

int
jac_stiff (double t, const double y[], double *dfdy, double dfdt[],
           void *params)
{
  gsl_matrix dfdy_mat;
  dfdy_mat.data = dfdy;
  dfdy_mat.size1 = 2;
  dfdy_mat.size2 = 2;
  dfdy_mat.tda = 2;
  dfdy_mat.block = 0;
  gsl_matrix_set (&dfdy_mat, 0, 0, 998.0);
  gsl_matrix_set (&dfdy_mat, 0, 1, 1998.0);
  gsl_matrix_set (&dfdy_mat, 1, 0, -999.0);
  gsl_matrix_set (&dfdy_mat, 1, 1, -1999.0);
  dfdt[0] = 0.0;
  dfdt[1] = 0.0;
  return GSL_SUCCESS;
}

gsl_odeiv_system rhs_func_stiff = {
  rhs_stiff,
  jac_stiff,
  2,
  0
};


void
test_stepper_linear (const gsl_odeiv_step_type * T, double h,
                     double base_prec)
{
  int s = 0;
  double y[2];
  double yerr[2];
  double t;
  double del;
  double delmax = 0.0;
  int count = 0;

  gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);

  y[0] = 1.0;
  y[1] = 0.0;

  for (t = 0.0; t < 4.0; t += h)
    {
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_lin);
      del = fabs ((y[1] - (t + h)) / y[1]);
      delmax = GSL_MAX_DBL (del, delmax);
      if (del > (count + 1.0) * base_prec)
        {
          printf ("  LINEAR(%20.17g)  %20.17g  %20.17g  %8.4g\n", t + h, y[1],
                  t + h, del);
          s++;
        }
      count++;
    }

  gsl_test (s, "%s, linear [0,4], max relative error = %g",
            gsl_odeiv_step_name (stepper), delmax);

  gsl_odeiv_step_free (stepper);
}


void
test_stepper_sin (const gsl_odeiv_step_type * T, double h, double base_prec)
{
  int s = 0;
  double y[2];
  double yerr[2];
  double t;
  double del;
  double delmax = 0.0;
  int count = 0;

  gsl_odeiv_step *stepper = gsl_odeiv_step_alloc (T, 2);

  y[0] = 1.0;
  y[1] = 0.0;

  for (t = 0.0; t < M_PI; t += h)
    {
      int stat;
      double sin_th = sin (t + h);
      gsl_odeiv_step_apply (stepper, t, h, y, yerr, 0, 0, &rhs_func_sin);
      del = fabs ((y[1] - sin_th) / sin_th);
      delmax = GSL_MAX_DBL (del, delmax);
      {
        if (t < 0.5 * M_PI)
          {
            stat = (del > (count + 1.0) * base_prec);
          }
        else if (t < 0.7 * M_PI)
          {
            stat = (del > 1.0e+04 * base_prec);
          }
        else if (t < 0.9 * M_PI)
          {
            stat = (del > 1.0e+06 * base_prec);
          }
        else
          {
            stat = (del > 1.0e+09 * base_prec);
          }
        if (stat != 0)
          {
            printf ("  SIN(%22.18g)  %22.18g  %22.18g  %10.6g\n", t + h, y[1],
                    sin_th, del);
          }
        s += stat;
      }
      count++;
    }
  if (delmax > 1.0e+09 * base_prec)
    {
      s++;
      printf ("  SIN(0 .. M_PI)  delmax = %g\n", delmax);
    }