アダムス法を用いた常微分方程式

ホーム > 製品 > nAG数値計算ライブラリ > サンプルソースコード集 > アダムス法を用いた常微分方程式 (C言語/C++)

Keyword: アダムス法, 常微分方程式

概要

本サンプルはアダムス法を用いた常微分方程式を求めるC言語によるサンプルプログラムです。本サンプルは以下に示される常微分方程式について４つのケースの計算をし、出力します。

アダムス法を用いた常微分方程式のデータ　

※本サンプルはnAG Cライブラリに含まれる関数 nag_ode_ivp_adams_gen() のExampleコードです。本サンプル及び関数の詳細情報は nag_ode_ivp_adams_gen のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで

出力結果

（本関数の詳細はnag_ode_ivp_adams_gen のマニュアルページを参照）

この出力例をダウンロード

nag_ode_ivp_adams_gen (d02cjc) Example Program Results

Case 1: intermediate output, root-finding

  Calculation with tol =    1.0e-04

     X         Y(1)         Y(2)         Y(3)
    0.00      0.50000      0.50000      0.62832
(User-supplied callback fcn, first invocation.)
(User-supplied callback g, first invocation.)
    2.00      1.54931      0.40548      0.30662
    4.00      1.74229      0.37433     -0.12890
    6.00      1.00554      0.41731     -0.55068

  Root of Y(1) = 0.0 at   7.288

  Solution is  -0.00000   0.47486  -0.76011

  Calculation with tol =    1.0e-05

     X         Y(1)         Y(2)         Y(3)
    0.00      0.50000      0.50000      0.62832
    2.00      1.54933      0.40548      0.30662
    4.00      1.74232      0.37433     -0.12891
    6.00      1.00552      0.41731     -0.55069

  Root of Y(1) = 0.0 at   7.288

  Solution is  -0.00000   0.47486  -0.76010


Case 2: no intermediate output, root-finding

  Calculation with tol =    1.0e-04

  Root of Y(1) = 0.0 at   7.288

  Solution is  -0.00000   0.47486  -0.76011

  Calculation with tol =    1.0e-05

  Root of Y(1) = 0.0 at   7.288

  Solution is  -0.00000   0.47486  -0.76010


Case 3: intermediate output, no root-finding

  Calculation with tol =    1.0e-04

     X         Y(1)         Y(2)         Y(3)
    0.00      0.50000      0.50000      0.62832
    2.00      1.54931      0.40548      0.30662
    4.00      1.74229      0.37433     -0.12890
    6.00      1.00554      0.41731     -0.55068
    8.00     -0.74589      0.51299     -0.85371
   10.00     -3.62813      0.63325     -1.05152

  Calculation with tol =    1.0e-05

     X         Y(1)         Y(2)         Y(3)
    0.00      0.50000      0.50000      0.62832
    2.00      1.54933      0.40548      0.30662
    4.00      1.74232      0.37433     -0.12891
    6.00      1.00552      0.41731     -0.55069
    8.00     -0.74601      0.51299     -0.85372
   10.00     -3.62829      0.63326     -1.05153


Case 4: no intermediate output, no root-finding ( integrate to xend)

  Calculation with tol =    1.0e-04

     X         Y(1)         Y(2)         Y(3)
    0.00      0.50000      0.50000      0.62832
   10.00     -3.62813      0.63325     -1.05152

  Calculation with tol =    1.0e-05

     X         Y(1)         Y(2)         Y(3)
    0.00      0.50000      0.50000      0.62832
   10.00     -3.62829      0.63326     -1.05153

３～２７行目にはケース１の結果が出力されています。ケース１ではy=0.0となるデータ点が見つかるまで 2.0の間隔でx=10.0まで計算を行い、中間結果を出力しています。
８～１１行目にxの値と許容値 1.0e-004 で計算したyの値が出力されています。
１３行目にy=0.0となるデータ点が出力されています。
１５行目に常微分方程式の解が出力されています。
２０～２３行目にxの値と許容値 1.0e-005 で計算したyの値が出力されています。
２５行目にy=0.0となるデータ点が出力されています。
２７行目に常微分方程式の解が出力されています。
３０～４２行目にはケース２の結果が出力されています。ケース２では中間結果の出力を行わず、解が見つかったら終了します。。
４５～６５行目にケース３の結果が出力されています。ケース３では中間結果を出力し、x=10.0まで計算を行っています。
６８～８０行目にケース４の結果が出力されています。ケース４では中間結果の出力を行わず、x=10.0まで計算を行っています。

ソースコード

（本関数の詳細はnag_ode_ivp_adams_gen のマニュアルページを参照）

※本サンプルソースコードはnAG数値計算ライブラリ（Windows, Linux, MAC等に対応）の関数を呼び出します。
サンプルのコンパイル及び実行方法

このソースコードをダウンロード

/* nag_ode_ivp_adams_gen (d02cjc) Example Program.
 *
 * CLL6I261D/CLL6I261DL Version.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 *
 */

#include <nag.h>
#include <math.h>
#include <stdio.h>
#include <nag_stdlib.h>
#include <nagd02.h>
#include <nagx01.h>

#ifdef __cplusplus
extern "C"
{
#endif
  static void nAG_CALL out(Integer neq, double *xsol, const double y[],
                           Nag_User *comm);
  static void nAG_CALL fcn(Integer neq, double x, const double y[],
                           double f[], Nag_User *comm);
  static double nAG_CALL g(Integer neq, double x, const double y[],
                           Nag_User *comm);
#ifdef __cplusplus
}
#endif

struct user
{
  double xend, h;
  Integer k;
  Integer *use_comm;
};

int main(void)
{
  static Integer use_comm[2] = { 1, 1 };
  Integer exit_status = 0, i, j, neq;
  Nag_User comm;
  double pi, tol, x, y[3];
  struct user s;
  NagError fail;

  INIT_FAIL(fail);

  printf("nag_ode_ivp_adams_gen (d02cjc) Example Program Results\n");

  /* For communication with user-supplied functions
   * assign address of user defined structure
   * to Nag pointer.
   */
  s.use_comm = use_comm;
  comm.p = (Pointer) &s;

  neq = 3;
  s.xend = 10.0;
  /* nag_pi (x01aac).
   * pi
   */
  pi = nag_pi;
  printf("\nCase 1: intermediate output, root-finding\n");
  for (j = 4; j <= 5; ++j) {
    tol = pow(10.0, (double) (-j));
    printf("\n  Calculation with tol = %10.1e\n", tol);
    x = 0.0;
    y[0] = 0.5;
    y[1] = 0.5;
    y[2] = pi / 5.0;
    s.k = 4;
    s.h = (s.xend - x) / (double) (s.k + 1);
    printf("\n     X         Y(1)         Y(2)         Y(3)\n");

    /* nag_ode_ivp_adams_gen (d02cjc).
     * Ordinary differential equation solver using a
     * variable-order variable-step Adams method (Black Box)
     */
    nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out, g,
                          &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    printf("\n  Root of Y(1) = 0.0 at %7.3f\n", x);
    printf("\n  Solution is");
    for (i = 0; i < 3; ++i)
      printf("%10.5f", y[i]);
    printf("\n");
  }
  printf("\n\nCase 2: no intermediate output, root-finding\n");
  for (j = 4; j <= 5; ++j) {
    tol = pow(10.0, (double) (-j));
    printf("\n  Calculation with tol = %10.1e\n", tol);
    x = 0.0;
    y[0] = 0.5;
    y[1] = 0.5;
    y[2] = pi / 5.0;

    /* nag_ode_ivp_adams_gen (d02cjc), see above. */
    nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, NULLFN, g,
                          &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
    printf("\n  Root of Y(1) = 0.0 at %7.3f\n", x);
    printf("\n  Solution is");
    for (i = 0; i < 3; ++i)
      printf("%10.5f", y[i]);
    printf("\n");
  }
  printf("\n\nCase 3: intermediate output, no root-finding\n");
  for (j = 4; j <= 5; ++j) {
    tol = pow(10.0, (double) (-j));
    printf("\n  Calculation with tol = %10.1e\n", tol);
    x = 0.0;
    y[0] = 0.5;
    y[1] = 0.5;
    y[2] = pi / 5.0;
    s.k = 4;
    s.h = (s.xend - x) / (double) (s.k + 1);
    printf("\n     X         Y(1)         Y(2)         Y(3)\n");

    /* nag_ode_ivp_adams_gen (d02cjc), see above. */
    nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, out,
                          NULLDFN, &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }
  }

  printf("\n\nCase 4: no intermediate output, no root-finding");
  printf(" ( integrate to xend)\n");
  for (j = 4; j <= 5; ++j) {
    tol = pow(10.0, (double) (-j));
    printf("\n  Calculation with tol = %10.1e\n", tol);
    x = 0.0;
    y[0] = 0.5;
    y[1] = 0.5;
    y[2] = pi / 5.0;
    printf("\n     X         Y(1)         Y(2)         Y(3)\n");
    printf("%8.2f", x);
    for (i = 0; i < 3; ++i)
      printf("%13.5f", y[i]);
    printf("\n");

    /* nag_ode_ivp_adams_gen (d02cjc), see above. */
    nag_ode_ivp_adams_gen(neq, fcn, &x, y, s.xend, tol, Nag_Mixed, NULLFN,
                          NULLDFN, &comm, &fail);
    if (fail.code != NE_NOERROR) {
      printf("Error from nag_ode_ivp_adams_gen (d02cjc).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    printf("%8.2f", x);
    for (i = 0; i < 3; ++i)
      printf("%13.5f", y[i]);
    printf("\n");
  }
END:
  return exit_status;
}

static void nAG_CALL out(Integer neq, double *xsol, const double y[],
                         Nag_User *comm)
{
  Integer i;
  struct user *s = (struct user *) comm->p;

  printf("%8.2f", *xsol);
  for (i = 0; i < 3; ++i) {
    printf("%13.5f", y[i]);
  }
  printf("\n");
  *xsol = s->xend - (double) s->k * s->h;
  s->k--;
}

static void nAG_CALL fcn(Integer neq, double x, const double y[], double f[],
                         Nag_User *comm)
{
  double pwr;
  struct user *s = (struct user *) comm->p;

  if (s->use_comm[0]) {
    printf("(User-supplied callback fcn, first invocation.)\n");
    s->use_comm[0] = 0;
  }

  f[0] = tan(y[2]);
  f[1] = -0.032 * tan(y[2]) / y[1] - 0.02 * y[1] / cos(y[2]);

  pwr = y[1];
  f[2] = -0.032 / (pwr * pwr);
}

static double nAG_CALL g(Integer neq, double x, const double y[],
                         Nag_User *comm)
{
  struct user *s = (struct user *) comm->p;

  if (s->use_comm[1]) {
    printf("(User-supplied callback g, first invocation.)\n");
    s->use_comm[1] = 0;
  }

  return y[0];
}

C言語によるサンプルソースコード : 使用関数名：nag_ode_ivp_adams_gen (d02cjc)

概要

出力結果

ソースコード