製品 & サービス 
サポート 
ダウンロード 
リソース 
ご購入 
会社案内 
お問い合わせ 

ケラーのボックス型スキームを用いた離散化による１階偏微分方程式

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

ホーム > 製品 > nAG数値計算ライブラリ > サンプルソースコード集 > ケラーのボックス型スキームを用いた離散化による１階偏微分方程式 (C言語/C++)

Keyword: ケラー, ボックス型スキーム, １階偏微分方程式

概要

本サンプルはケラーのボックス型スキームを用いた離散化により１階偏微分方程式を解くC言語によるサンプルプログラムです。 本サンプルは以下に示される１階偏微分方程式を解いて結果を出力します。

１階偏微分方程式のデータ　

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

出力結果

（本関数の詳細はnag_pde_parab_1d_keller のマニュアルページを参照） 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
この出力例をダウンロード 
nag_pde_parab_1d_keller (d03pec) Example Program Results

  Accuracy requirement =   1.000e-06 Number of points =  41

 x            0.1000    0.3000    0.5000    0.7000    0.9000

(User-supplied callback bndary, first invocation.)
(User-supplied callback pdedef, first invocation.)
 t =  0.20
 Approx u1    0.7845    1.0010    1.2733    1.6115    2.0281
 Exact  u1    0.7845    1.0010    1.2733    1.6115    2.0281
 Approx u2   -0.8352   -0.8159   -0.8367   -0.9128   -1.0609
 Exact  u2   -0.8353   -0.8160   -0.8367   -0.9129   -1.0609

 t =  0.40
 Approx u1    0.6481    0.8533    1.1212    1.4627    1.8903
 Exact  u1    0.6481    0.8533    1.1212    1.4627    1.8903
 Approx u2   -1.5216   -1.6767   -1.8934   -2.1917   -2.5944
 Exact  u2   -1.5217   -1.6767   -1.8935   -2.1917   -2.5945

 t =  0.60
 Approx u1    0.6892    0.8961    1.1747    1.5374    1.9989
 Exact  u1    0.6892    0.8962    1.1747    1.5374    1.9989
 Approx u2   -2.0047   -2.3434   -2.7677   -3.3002   -3.9680
 Exact  u2   -2.0048   -2.3436   -2.7678   -3.3003   -3.9680

 t =  0.80
 Approx u1    0.8977    1.1247    1.4320    1.8349    2.3514
 Exact  u1    0.8977    1.1247    1.4320    1.8349    2.3512
 Approx u2   -2.3403   -2.8675   -3.5110   -4.2960   -5.2536
 Exact  u2   -2.3405   -2.8677   -3.5111   -4.2961   -5.2537

 t =  1.00
 Approx u1    1.2470    1.5206    1.8828    2.3528    2.9519
 Exact  u1    1.2470    1.5205    1.8829    2.3528    2.9518
 Approx u2   -2.6229   -3.3338   -4.1998   -5.2505   -6.5218
 Exact  u2   -2.6232   -3.3340   -4.2001   -5.2507   -6.5219

 Number of integration steps in time =    149
 Number of function evaluations =    399
 Number of Jacobian evaluations =    13
 Number of iterations =    323

  • ３行目には誤差推定を制御するための正数とメッシュ点の数が出力されています。
  • ５行目には空間方向のメッシュ点が出力されています。
  • ７～１１行目にはt=0.20の場合のU1(x,t)の近似解と厳密解、U2(x,t)の近似解と厳密解が出力されています。
  • １３～１７行目にはt=0.40の場合のU1(x,t)の近似解と厳密解、U2(x,t)の近似解と厳密解が出力されています。
  • １９～２３行目にはt=0.60の場合のU1(x,t)の近似解と厳密解、U2(x,t)の近似解と厳密解が出力されています。
  • ２５～２９行目にはt=0.80の場合のU1(x,t)の近似解と厳密解、U2(x,t)の近似解と厳密解が出力されています。
  • ３１～３５行目にはt=1.00の場合のU1(x,t)の近似解と厳密解、U2(x,t)の近似解と厳密解が出力されています。
  • ３７行目には積分のステップの数が出力されています。
  • ３８行目には関数評価の数が出力されています。
  • ３９行目にはヤコビアン評価の数が出力されています。
  • ４０行目には反復数が出力されています。

ソースコード

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

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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
このソースコードをダウンロード 
/* nag_pde_parab_1d_keller (d03pec) Example Program.
 *
 * CLL6I261D/CLL6I261DL Version.
 *
 * Copyright 2017 Numerical Algorithms Group.
 *
 * Mark 26.1, 2017.
 */

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

#ifdef __cplusplus
extern "C"
{
#endif
  static void nAG_CALL pdedef(Integer, double, double, const double[],
                              const double[], const double[], double[],
                              Integer *, Nag_Comm *);
  static void nAG_CALL bndary(Integer, double, Integer, Integer,
                              const double[], const double[], double[],
                              Integer *, Nag_Comm *);
  static void nAG_CALL exact(double, Integer, Integer, double *, double *);
  static void nAG_CALL uinit(Integer, Integer, double *, double *);
#ifdef __cplusplus
}
#endif

#define U(I, J)  u[npde*((J) -1)+(I) -1]
#define EU(I, J) eu[npde*((J) -1)+(I) -1]

int main(void)
{
  const Integer npde = 2, npts = 41, nleft = 1, neqn = npde * npts;
  const Integer lisave = neqn + 24, nwkres =
         npde * (npts + 21 + 3 * npde) + 7 * npts + 4;
  const Integer lrsave =
         11 * neqn + (4 * npde + nleft + 2) * neqn + 50 + nwkres;
  static double ruser[2] = { -1.0, -1.0 };
  Integer exit_status = 0, i, ind, it, itask, itrace;
  double acc, tout, ts;
  double *eu = 0, *rsave = 0, *u = 0, *x = 0;
  Integer *isave = 0;
  NagError fail;
  Nag_Comm comm;
  Nag_D03_Save saved;

  INIT_FAIL(fail);

  printf("nag_pde_parab_1d_keller (d03pec) Example Program Results\n\n");

  /* For communication with user-supplied functions: */
  comm.user = ruser;

  /* Allocate memory */

  if (!(eu = nAG_ALLOC(npde * npts, double)) ||
      !(rsave = nAG_ALLOC(lrsave, double)) ||
      !(u = nAG_ALLOC(npde * npts, double)) ||
      !(x = nAG_ALLOC(npts, double)) || !(isave = nAG_ALLOC(lisave, Integer)))
  {
    printf("Allocation failure\n");
    exit_status = 1;
    goto END;
  }

  itrace = 0;
  acc = 1e-6;

  printf("  Accuracy requirement =%12.3e", acc);
  printf(" Number of points = %3ld\n\n", npts);

  /* Set spatial-mesh points */

  for (i = 0; i < npts; ++i)
    x[i] = i / (npts - 1.0);

  printf(" x        ");
  printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n\n",
         x[4], x[12], x[20], x[28], x[36]);

  ind = 0;
  itask = 1;

  uinit(npde, npts, x, u);

  /* Loop over output value of t */

  ts = 0.0;
  for (it = 0; it < 5; ++it) {
    tout = 0.2 * (it + 1);
    /* nag_pde_parab_1d_keller (d03pec).
     * General system of first-order PDEs, method of lines,
     * Keller box discretization, one space variable
     */
    nag_pde_parab_1d_keller(npde, &ts, tout, pdedef, bndary, u, npts, x,
                            nleft, acc, rsave, lrsave, isave, lisave, itask,
                            itrace, 0, &ind, &comm, &saved, &fail);

    if (fail.code != NE_NOERROR) {
      printf("Error from nag_pde_parab_1d_keller (d03pec).\n%s\n",
             fail.message);
      exit_status = 1;
      goto END;
    }

    /* Check against the exact solution */

    exact(tout, npde, npts, x, eu);

    printf(" t = %5.2f\n", ts);
    printf(" Approx u1");
    printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n",
           U(1, 5), U(1, 13), U(1, 21), U(1, 29), U(1, 37));

    printf(" Exact  u1");
    printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n",
           EU(1, 5), EU(1, 13), EU(1, 21), EU(1, 29), EU(1, 37));

    printf(" Approx u2");
    printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n",
           U(2, 5), U(2, 13), U(2, 21), U(2, 29), U(2, 37));

    printf(" Exact  u2");
    printf("%10.4f%10.4f%10.4f%10.4f%10.4f\n\n",
           EU(2, 5), EU(2, 13), EU(2, 21), EU(2, 29), EU(2, 37));
  }
  printf(" Number of integration steps in time = %6ld\n", isave[0]);
  printf(" Number of function evaluations = %6ld\n", isave[1]);
  printf(" Number of Jacobian evaluations =%6ld\n", isave[2]);
  printf(" Number of iterations = %6ld\n\n", isave[4]);

END:
  nAG_FREE(eu);
  nAG_FREE(rsave);
  nAG_FREE(u);
  nAG_FREE(x);
  nAG_FREE(isave);

  return exit_status;
}

static void nAG_CALL pdedef(Integer npde, double t, double x,
                            const double u[], const double udot[],
                            const double dudx[], double res[], Integer *ires,
                            Nag_Comm *comm)
{
  if (comm->user[0] == -1.0) {
    printf("(User-supplied callback pdedef, first invocation.)\n");
    comm->user[0] = 0.0;
  }
  if (*ires == -1) {
    res[0] = udot[0];
    res[1] = udot[1];
  }
  else {
    res[0] = udot[0] + dudx[0] + dudx[1];
    res[1] = udot[1] + 4.0 * dudx[0] + dudx[1];
  }
  return;
}

static void nAG_CALL bndary(Integer npde, double t, Integer ibnd,
                            Integer nobc, const double u[],
                            const double udot[], double res[], Integer *ires,
                            Nag_Comm *comm)
{
  if (comm->user[1] == -1.0) {
    printf("(User-supplied callback bndary, first invocation.)\n");
    comm->user[1] = 0.0;
  }
  if (ibnd == 0) {
    if (*ires == -1) {
      res[0] = 0.0;
    }
    else {
      res[0] = u[0] - 0.5 * (exp(t) + exp(-3.0 * t))
             - 0.25 * (sin(-3.0 * t) - sin(t));
    }
  }
  else {
    if (*ires == -1) {
      res[0] = 0.0;
    }
    else {
      res[0] = u[1] - exp(1.0 - 3.0 * t) + exp(t + 1.0)
             - 0.5 * (sin(1.0 - 3.0 * t) + sin(t + 1.0));
    }
  }
  return;
}

static void nAG_CALL uinit(Integer npde, Integer npts, double *x, double *u)
{
  /* Routine for PDE initial values */

  Integer i;

  for (i = 1; i <= npts; ++i) {
    U(1, i) = exp(x[i - 1]);
    U(2, i) = sin(x[i - 1]);
  }
  return;
}

static void nAG_CALL exact(double t, Integer npde, Integer npts, double *x,
                           double *u)
{
  /* Exact solution (for comparison purposes) */

  Integer i;

  for (i = 1; i <= npts; ++i) {
    U(1, i) = 0.5 * (exp(x[i - 1] + t) + exp(x[i - 1] - 3.0 * t)) +
           0.25 * (sin(x[i - 1] - 3.0 * t) - sin(x[i - 1] + t));
    U(2, i) = exp(x[i - 1] - 3.0 * t) - exp(x[i - 1] + t) +
           0.5 * (sin(x[i - 1] - 3.0 * t) + sin(x[i - 1] + t));
  }
  return;
}


関連情報
© 日本ニューメリカルアルゴリズムズグループ株式会社 2025
Privacy Policy  /  Trademarks