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