Keyword: ケラー, ボックス型スキーム, １階偏微分方程式
概要
本サンプルはケラーのボックス型スキームを用いた離散化により１階偏微分方程式を解くFortranによるサンプルプログラムです。 本サンプルは以下に示される１階偏微分方程式を解いて結果を出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン d03pef() のExampleコードです。本サンプル及びルーチンの詳細情報は d03pef のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はd03pef のマニュアルページを参照）- １行目はタイトル行で読み飛ばされます。
- ２行目はメッシュ点の数（npts=41）を指定しています。
- ３行目は局所誤差を制御するための正数（acc=0.1E-5）を指定しています。
- ４行目はトレース情報のレベル（itrace=0:警告メッセージのみ出力）を指定しています。
- ５行目は独立変数tの初期値（ts=0.0）と積分が実行されるtの最終値（tout=0.0）を指定しています。
出力結果
（本ルーチンの詳細はd03pef のマニュアルページを参照）| この出力例をダウンロード |
D03PEF Example Program Results Accuracy requirement = 0.100E-05 Number of points = 41 X 0.1000 0.3000 0.5000 0.7000 0.9000 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)の近似解と厳密解が出力されています。
- ３８行目には積分のステップの数が出力されています。
- ３９行目には関数評価の数が出力されています。
- ４０行目にはヤコビアン評価の数が出力されています。
- ４１行目には反復数が出力されています。
ソースコード
（本ルーチンの詳細はd03pef のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「nAG Fortranライブラリ」のルーチンを呼び出します。
サンプルのコンパイル及び実行方法
| このソースコードをダウンロード |
! D03PEF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
MODULE d03pefe_mod
! D03PEF Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
USE nag_library, ONLY : nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
REAL (KIND=nag_wp), PARAMETER :: half = 0.5_nag_wp
REAL (KIND=nag_wp), PARAMETER :: one = 1.0_nag_wp
INTEGER, PARAMETER :: nin = 5, nleft = 1, nout = 6, &
npde = 2
CONTAINS
SUBROUTINE uinit(npde,npts,x,u)
! Routine for PDE initial values
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
INTEGER, INTENT (IN) :: npde, npts
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: u(npde,npts)
REAL (KIND=nag_wp), INTENT (IN) :: x(npts)
! .. Local Scalars ..
INTEGER :: i
! .. Intrinsic Functions ..
INTRINSIC exp, sin
! .. Executable Statements ..
DO i = 1, npts
u(1,i) = exp(x(i))
u(2,i) = sin(x(i))
END DO
RETURN
END SUBROUTINE uinit
SUBROUTINE pdedef(npde,t,x,u,ut,ux,res,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t, x
INTEGER, INTENT (INOUT) :: ires
INTEGER, INTENT (IN) :: npde
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: res(npde)
REAL (KIND=nag_wp), INTENT (IN) :: u(npde), ut(npde), ux(npde)
! .. Executable Statements ..
IF (ires==-1) THEN
res(1) = ut(1)
res(2) = ut(2)
ELSE
res(1) = ut(1) + ux(1) + ux(2)
res(2) = ut(2) + 4.0_nag_wp*ux(1) + ux(2)
END IF
RETURN
END SUBROUTINE pdedef
SUBROUTINE bndary(npde,t,ibnd,nobc,u,ut,res,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t
INTEGER, INTENT (IN) :: ibnd, nobc, npde
INTEGER, INTENT (INOUT) :: ires
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: res(nobc)
REAL (KIND=nag_wp), INTENT (IN) :: u(npde), ut(npde)
! .. Local Scalars ..
REAL (KIND=nag_wp) :: t1, t3
! .. Intrinsic Functions ..
INTRINSIC exp, sin
! .. Executable Statements ..
IF (ires==-1) THEN
res(1) = 0.0_nag_wp
ELSE IF (ibnd==0) THEN
t3 = -3.0_nag_wp*t
t1 = t
res(1) = u(1) - half*((exp(t3)+exp(t1))+half*(sin(t3)-sin(t1)))
ELSE
t3 = one - 3.0_nag_wp*t
t1 = one + t
res(1) = u(2) - ((exp(t3)-exp(t1))+half*(sin(t3)+sin(t1)))
END IF
RETURN
END SUBROUTINE bndary
SUBROUTINE exact(t,npde,npts,x,u)
! Exact solution (for comparison purposes)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t
INTEGER, INTENT (IN) :: npde, npts
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: u(npde,npts)
REAL (KIND=nag_wp), INTENT (IN) :: x(npts)
! .. Local Scalars ..
REAL (KIND=nag_wp) :: xt, xt3
INTEGER :: i
! .. Intrinsic Functions ..
INTRINSIC exp, sin
! .. Executable Statements ..
DO i = 1, npts
xt3 = x(i) - 3.0_nag_wp*t
xt = x(i) + t
u(1,i) = half*((exp(xt3)+exp(xt))+half*(sin(xt3)-sin(xt)))
u(2,i) = (exp(xt3)-exp(xt)) + half*(sin(xt3)+sin(xt))
END DO
RETURN
END SUBROUTINE exact
END MODULE d03pefe_mod
PROGRAM d03pefe
! D03PEF Example Main Program
! .. Use Statements ..
USE nag_library, ONLY : d03pef, nag_wp
USE d03pefe_mod, ONLY : bndary, exact, nin, nleft, nout, npde, pdedef, &
uinit
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Local Scalars ..
REAL (KIND=nag_wp) :: acc, tout, ts
INTEGER :: i, ifail, ind, it, itask, &
itrace, lisave, lrsave, neqn, &
npts, nwkres
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: eu(:,:), rsave(:), u(:,:), x(:)
INTEGER, ALLOCATABLE :: isave(:)
! .. Intrinsic Functions ..
INTRINSIC real
! .. Executable Statements ..
WRITE (nout,*) 'D03PEF Example Program Results'
! Skip heading in data file
READ (nin,*)
READ (nin,*) npts
nwkres = npde*(npts+21+3*npde) + 7*npts + 4
neqn = npde*npts
lisave = neqn + 24
lrsave = 11*neqn + (4*npde+nleft+2)*neqn + 50 + nwkres
ALLOCATE (eu(npde,npts),rsave(lrsave),u(npde,npts),x(npts), &
isave(lisave))
READ (nin,*) acc
READ (nin,*) itrace
! Set spatial-mesh points
DO i = 1, npts
x(i) = real(i-1,kind=nag_wp)/real(npts-1,kind=nag_wp)
END DO
ind = 0
itask = 1
CALL uinit(npde,npts,x,u)
! Loop over output value of t
READ (nin,*) ts, tout
DO it = 1, 5
tout = 0.2_nag_wp*real(it,kind=nag_wp)
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
CALL d03pef(npde,ts,tout,pdedef,bndary,u,npts,x,nleft,acc,rsave, &
lrsave,isave,lisave,itask,itrace,ind,ifail)
IF (it==1) THEN
WRITE (nout,99997) acc, npts
WRITE (nout,99999) x(5), x(13), x(21), x(29), x(37)
END IF
! Check against the exact solution
CALL exact(tout,npde,npts,x,eu)
WRITE (nout,99998) ts
WRITE (nout,99995) u(1,5:37:8)
WRITE (nout,99994) eu(1,5:37:8)
WRITE (nout,99993) u(2,5:37:8)
WRITE (nout,99992) eu(2,5:37:8)
END DO
WRITE (nout,99996) isave(1), isave(2), isave(3), isave(5)
99999 FORMAT (' X ',5F10.4/)
99998 FORMAT (' T = ',F5.2)
99997 FORMAT (//' Accuracy requirement =',E10.3,' Number of points = ',I3/)
99996 FORMAT (' Number of integration steps in time = ',I6/' Number o', &
'f function evaluations = ',I6/' Number of Jacobian eval', &
'uations =',I6/' Number of iterations = ',I6)
99995 FORMAT (' Approx U1',5F10.4)
99994 FORMAT (' Exact U1',5F10.4)
99993 FORMAT (' Approx U2',5F10.4)
99992 FORMAT (' Exact U2',5F10.4/)
END PROGRAM d03pefe