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

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

Fortranによるサンプルソースコード : 使用ルーチン名：d02cjf

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

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

概要

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

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

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

入力データ

（本ルーチンの詳細はd02cjf のマニュアルページを参照） 
1
2
3
4
5
このデータをダウンロード 
D02CJF Example Program Data
   0.0                                    : xinit
  10.0                                    : xend
   0.5  0.5  6.28318530717958647692E-1    : yinit
   4                                      : kinit

  • １行目はタイトル行で読み飛ばされます。
  • ２行目は独立変数xの初期値(xinit)を指定しています。　
  • ３行目は独立変数xの最終値(xend)を指定しています。
  • ４行目は解の初期値(yinit)を指定しています。
  • ５行目は独立変数xの刻み幅の値の数(kinit)を指定しています。

出力結果

（本ルーチンの詳細はd02cjf のマニュアルページを参照） 
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
この出力例をダウンロード 
 D02CJF Example Program Results

 Case 1: intermediate output, root-finding

  Calculation with TOL = 0.1E-03
      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
   Root of Y(1) = 0.0 at  7.288
   Solution is      0.00000      0.47486     -0.76011

  Calculation with TOL = 0.1E-04
      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 = 0.1E-03
   Root of Y(1) = 0.0 at  7.288
   Solution is      0.00000      0.47486     -0.76011

  Calculation with TOL = 0.1E-04
   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 = 0.1E-03
      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 = 0.1E-04
      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 = 0.1E-03
      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 = 0.1E-04
      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の値と許容値 0.1e-003 で計算したyの値が出力されています。
  • １１行目にy=0.0となるデータ点が出力されています。
  • １２行目に常微分方程式の解が出力されています。
  • １４～１９行目にxの値と許容値 0.1e-004 で計算したyの値が出力されています。
  • ２０行目にy=0.0となるデータ点が出力されています。
  • ２１行目に常微分方程式の解が出力されています。
  • ２４～３２行目にはケース２の結果が出力されています。ケース２では中間結果の出力を行わず、解が見つかったら終了します。。
  • ３５～５３行目にケース３の結果が出力されています。ケース３では中間結果を出力し、x=10.0まで計算を行っています。
  • ５６～６６行目にケース４の結果が出力されています。ケース４では中間結果の出力を行わず、x=10.0まで計算を行っています。

ソースコード

（本ルーチンの詳細はd02cjf のマニュアルページを参照）

※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法

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
このソースコードをダウンロード 
!   D02CJF Example Program Text
!   Mark 23 Release. nAG Copyright 2011.

    MODULE d02cjfe_mod

!      Data for D02CJF example program

!      .. Use Statements ..
       USE nag_library, ONLY : nag_wp
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Parameters ..
       INTEGER, PARAMETER                  :: n = 3, nin = 5, nout = 6
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: h, xend
       INTEGER, SAVE                       :: k
!   n: number of differential equations
    CONTAINS
       SUBROUTINE output(xsol,y)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (INOUT)  :: xsol
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: y(*)
!         .. Local Scalars ..
          INTEGER                             :: j
!         .. Intrinsic Functions ..
          INTRINSIC                              real
!         .. Executable Statements ..
          WRITE (nout,99999) xsol, (y(j),j=1,n)
          xsol = xend - real(k,kind=nag_wp)*h
          k = k - 1
          RETURN

99999     FORMAT (1X,F8.2,3F13.5)
       END SUBROUTINE output
       SUBROUTINE fcn(x,y,f)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Parameters ..
          REAL (KIND=nag_wp), PARAMETER       :: alpha = -0.032E0_nag_wp
          REAL (KIND=nag_wp), PARAMETER       :: beta = -0.02E0_nag_wp
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (OUT)    :: f(*)
          REAL (KIND=nag_wp), INTENT (IN)     :: y(*)
!         .. Intrinsic Functions ..
          INTRINSIC                              cos, tan
!         .. Executable Statements ..
          f(1) = tan(y(3))
          f(2) = alpha*tan(y(3))/y(2) + beta*y(2)/cos(y(3))
          f(3) = alpha/y(2)**2
          RETURN
       END SUBROUTINE fcn
       FUNCTION g(x,y)

!         .. Implicit None Statement ..
          IMPLICIT NONE
!         .. Function Return Value ..
          REAL (KIND=nag_wp)                  :: g
!         .. Scalar Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: x
!         .. Array Arguments ..
          REAL (KIND=nag_wp), INTENT (IN)     :: y(*)
!         .. Executable Statements ..
          g = y(1)
          RETURN
       END FUNCTION g
    END MODULE d02cjfe_mod
    PROGRAM d02cjfe

!      D02CJF Example Main Program

!      .. Use Statements ..
       USE nag_library, ONLY : d02cjf, d02cjw, d02cjx, nag_wp
       USE d02cjfe_mod, ONLY : fcn, g, h, k, n, nin, nout, output, xend
!      .. Implicit None Statement ..
       IMPLICIT NONE
!      .. Local Scalars ..
       REAL (KIND=nag_wp)                  :: tol, x, xinit
       INTEGER                             :: i, icase, ifail, iw, j, kinit
!      .. Local Arrays ..
       REAL (KIND=nag_wp), ALLOCATABLE     :: w(:), y(:), yinit(:)
!      .. Intrinsic Functions ..
       INTRINSIC                              real
!      .. Executable Statements ..
       WRITE (nout,*) 'D02CJF Example Program Results'
       iw = 21*n + 28
       ALLOCATE (w(iw),y(n),yinit(n))
!      Skip heading in data file
       READ (nin,*)
!      xinit: initial x value, xend: final x value.
       READ (nin,*) xinit
       READ (nin,*) xend
       READ (nin,*) yinit(1:n)
       READ (nin,*) kinit
       DO icase = 1, 4
          WRITE (nout,*)
          SELECT CASE (icase)
          CASE (1)
             WRITE (nout,99995) icase, 'intermediate output, root-finding'
          CASE (2)
             WRITE (nout,99995) icase, 'no intermediate output, root-finding'
          CASE (3)
             WRITE (nout,99995) icase, 'intermediate output, no root-finding'
          CASE (4)
             WRITE (nout,99995) icase, 'no intermediate output, &
                &no root-finding ( integrate to XEND)'
          END SELECT
          DO j = 4, 5
             tol = 10.0E0_nag_wp**(-j)
             WRITE (nout,*)
             WRITE (nout,99999) ' Calculation with TOL =', tol
             x = xinit
             y(1:n) = yinit(1:n)
             IF (icase/=2) THEN
                WRITE (nout,*) '     X         Y(1)         Y(2)         Y(3)'
                k = kinit
                h = (xend-x)/real(k+1,kind=nag_wp)
             END IF
!            ifail: behaviour on error exit   
!                   =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
             ifail = 0
             SELECT CASE (icase)
             CASE (1)
                CALL d02cjf(x,xend,n,y,fcn,tol,'Default',output,g,w,ifail)
                WRITE (nout,99998) '  Root of Y(1) = 0.0 at', x
                WRITE (nout,99997) '  Solution is', (y(i),i=1,n)
             CASE (2)
                CALL d02cjf(x,xend,n,y,fcn,tol,'Default',d02cjx,g,w,ifail)
                WRITE (nout,99998) '  Root of Y(1) = 0.0 at', x
                WRITE (nout,99997) '  Solution is', (y(i),i=1,n)
             CASE (3)
                CALL d02cjf(x,xend,n,y,fcn,tol,'Default',output,d02cjw,w, &
                   ifail)
             CASE (4)
                WRITE (nout,99996) x, (y(i),i=1,n)
                CALL d02cjf(x,xend,n,y,fcn,tol,'Default',d02cjx,d02cjw,w, &
                   ifail)
                WRITE (nout,99996) x, (y(i),i=1,n)
             END SELECT
          END DO
          IF (icase<4) THEN
             WRITE (nout,*)
          END IF
       END DO

99999  FORMAT (1X,A,E8.1)
99998  FORMAT (1X,A,F7.3)
99997  FORMAT (1X,A,3F13.5)
99996  FORMAT (1X,F8.2,3F13.5)
99995  FORMAT (1X,'Case ',I1,': ',A)
    END PROGRAM d02cjfe


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