Keyword: アダムス法, 常微分方程式
概要
本サンプルはアダムス法を用いた常微分方程式を求めるFortranによるサンプルプログラムです。 本サンプルは以下に示される常微分方程式について４つのケースの計算をし、出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン d02cjf() のExampleコードです。本サンプル及びルーチンの詳細情報は d02cjf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はd02cjf のマニュアルページを参照）| このデータをダウンロード |
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 のマニュアルページを参照）| この出力例をダウンロード |
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ライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
! 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