Keyword: チェビシェフ, 選点法, 放物型偏微分方程式
概要
本サンプルはチェビシェフC0選点法を用いた離散化により放物型偏微分方程式を解くFortranによるサンプルプログラムです。 本サンプルは以下に示される１組の２次の放物型偏微分方程式（PDE）で記述される４次の偏微分方程式（PDE）を解いて結果を出力します。
※本サンプルはnAG Fortranライブラリに含まれるルーチン d03pdf() のExampleコードです。本サンプル及びルーチンの詳細情報は d03pdf のマニュアルページをご参照ください。
ご相談やお問い合わせはこちらまで
入力データ
（本ルーチンの詳細はd03pdf のマニュアルページを参照）1 2 3 4 5 6
| このデータをダウンロード |
D03PDF Example Program Data 6 10 3 1 : intpts, nbkpts, npoly, itype -1.0 -0.6 -0.2 0.2 0.6 1.0 : xout 1.0E-4 : acc 0 0 : m, itrace 0.0 0.1E-4 : ts, tout
- １行目はタイトル行で読み飛ばされます。
- ２行目は補間数（intpts=6）、ブレークポイントの数（nbkpts=10）、偏微分方程式（PDE）の解の近似に使用されるチェビシェフ多項式の次数（npoly=3）、実行される補間（itype=1:補間点の解を計算）を指定しています。
- ３行目は空間補間点(xout)を指定しています。
- ４行目は局所誤差を制御するための正数(acc)を指定しています。
- ５行目は座標システム（m=0:デカルト座標）とトレース情報のレベル（itrace=0:警告メッセージのみ出力）を指定しています。
- ６行目は独立変数tの初期値(ts)と積分が実行されるtの最終値(tout)を指定しています。
出力結果
（本ルーチンの詳細はd03pdf のマニュアルページを参照）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
| この出力例をダウンロード |
D03PDF Example Program Results
Polynomial degree = 3 No. of elements = 9
Accuracy requirement = 0.100E-03 Number of points = 28
T / X -1.0000 -0.6000 -0.2000 0.2000 0.6000 1.0000
0.0001 U(1) 1.0000 0.8090 0.3090 -0.3090 -0.8090 -1.0000
U(2) -2.4850 -1.9957 -0.7623 0.7623 1.9957 2.4850
0.0010 U(1) 1.0000 0.8085 0.3088 -0.3088 -0.8085 -1.0000
U(2) -2.5583 -1.9913 -0.7606 0.7606 1.9913 2.5583
0.0100 U(1) 1.0000 0.8051 0.3068 -0.3068 -0.8051 -1.0000
U(2) -2.6962 -1.9481 -0.7439 0.7439 1.9481 2.6962
0.1000 U(1) 1.0000 0.7951 0.2985 -0.2985 -0.7951 -1.0000
U(2) -2.9022 -1.8339 -0.6338 0.6338 1.8339 2.9022
1.0000 U(1) 1.0000 0.7939 0.2972 -0.2972 -0.7939 -1.0000
U(2) -2.9233 -1.8247 -0.6120 0.6120 1.8247 2.9233
Number of integration steps in time 50
Number of residual evaluations of resulting ODE system 407
Number of Jacobian evaluations 18
Number of iterations of nonlinear solver 122
- ２行目にはチェビシェフ多項式の次数と要素の数が出力されています。
- ３行目には誤差推定を制御するための正数とメッシュ点の数が出力されています。
- ５行目にはブレークポイントが出力されています。
- ７～２０行目には独立変数の値とU1(x,t)とU2(x,t)の計算結果が出力されています。
- ２２行目には積分のステップの数が出力されています。
- ２３行目には結果として生じる常微分方程式（ODE）の残差評価の数が出力されています。
- ２４行目にはヤコビアン評価の数が出力されています。
- ２５行目には非線形ソルバの反復数が出力されています。
ソースコード
（本ルーチンの詳細はd03pdf のマニュアルページを参照）
※本サンプルソースコードは科学技術・統計計算ライブラリである「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 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
| このソースコードをダウンロード |
! D03PDF Example Program Text
! Mark 23 Release. nAG Copyright 2011.
MODULE d03pdfe_mod
! D03PDF Example Program Module:
! Parameters and User-defined Routines
! .. Use Statements ..
USE nag_library, ONLY : nag_wp
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Parameters ..
INTEGER, PARAMETER :: nin = 5, nout = 6, npde = 2
CONTAINS
SUBROUTINE uinit(npde,npts,x,u)
! .. Use Statements ..
USE nag_library, ONLY : x01aaf
! .. 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 ..
REAL (KIND=nag_wp) :: piby2
INTEGER :: i
! .. Intrinsic Functions ..
INTRINSIC sin
! .. Executable Statements ..
piby2 = 0.5_nag_wp*x01aaf(piby2)
DO i = 1, npts
u(1,i) = -sin(piby2*x(i))
u(2,i) = -piby2*piby2*u(1,i)
END DO
RETURN
END SUBROUTINE uinit
SUBROUTINE pdedef(npde,t,x,nptl,u,ux,p,q,r,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t
INTEGER, INTENT (INOUT) :: ires
INTEGER, INTENT (IN) :: npde, nptl
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: p(npde,npde,nptl), &
q(npde,nptl), r(npde,nptl)
REAL (KIND=nag_wp), INTENT (IN) :: u(npde,nptl), ux(npde,nptl), &
x(nptl)
! .. Local Scalars ..
INTEGER :: i
! .. Executable Statements ..
DO i = 1, nptl
q(1,i) = u(2,i)
q(2,i) = u(1,i)*ux(2,i) - ux(1,i)*u(2,i)
r(1,i) = ux(1,i)
r(2,i) = ux(2,i)
p(1,1,i) = 0.0_nag_wp
p(1,2,i) = 0.0_nag_wp
p(2,1,i) = 0.0_nag_wp
p(2,2,i) = 1.0_nag_wp
END DO
RETURN
END SUBROUTINE pdedef
SUBROUTINE bndary(npde,t,u,ux,ibnd,beta,gamma,ires)
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Scalar Arguments ..
REAL (KIND=nag_wp), INTENT (IN) :: t
INTEGER, INTENT (IN) :: ibnd, npde
INTEGER, INTENT (INOUT) :: ires
! .. Array Arguments ..
REAL (KIND=nag_wp), INTENT (OUT) :: beta(npde), gamma(npde)
REAL (KIND=nag_wp), INTENT (IN) :: u(npde), ux(npde)
! .. Executable Statements ..
IF (ibnd==0) THEN
beta(1) = 1.0_nag_wp
gamma(1) = 0.0_nag_wp
beta(2) = 0.0_nag_wp
gamma(2) = u(1) - 1.0_nag_wp
ELSE
beta(1) = 1.0E+0_nag_wp
gamma(1) = 0.0_nag_wp
beta(2) = 0.0_nag_wp
gamma(2) = u(1) + 1.0_nag_wp
END IF
RETURN
END SUBROUTINE bndary
END MODULE d03pdfe_mod
PROGRAM d03pdfe
! D03PDF Example Main Program
! .. Use Statements ..
USE nag_library, ONLY : d03pdf, d03pyf, nag_wp
USE d03pdfe_mod, ONLY : bndary, nin, nout, npde, pdedef, uinit
! .. Implicit None Statement ..
IMPLICIT NONE
! .. Local Scalars ..
REAL (KIND=nag_wp) :: acc, dx, tout, ts
INTEGER :: i, ifail, ind, intpts, it, &
itask, itrace, itype, lenode, &
lisave, lrsave, m, mu, nbkpts, &
nel, neqn, npl1, npoly, npts, &
nwkres
! .. Local Arrays ..
REAL (KIND=nag_wp), ALLOCATABLE :: rsave(:), u(:,:), uout(:,:,:), &
x(:), xbkpts(:), xout(:)
INTEGER, ALLOCATABLE :: isave(:)
! .. Intrinsic Functions ..
INTRINSIC real
! .. Executable Statements ..
WRITE (nout,*) 'D03PDF Example Program Results'
! Skip heading in data file
READ (nin,*)
READ (nin,*) intpts, nbkpts, npoly, itype
nel = nbkpts - 1
npts = nel*npoly + 1
mu = npde*(npoly+1) - 1
neqn = npde*npts
lisave = neqn + 24
npl1 = npoly + 1
nwkres = 3*npl1*npl1 + npl1*(npde*npde+6*npde+nbkpts+1) + 13*npde + 5
lenode = (3*mu+1)*neqn
lrsave = 11*neqn + 50 + nwkres + lenode
ALLOCATE (u(npde,npts),uout(npde,intpts,itype),rsave(lrsave),x(npts), &
xbkpts(nbkpts),xout(intpts),isave(lisave))
READ (nin,*) xout(1:intpts)
READ (nin,*) acc
READ (nin,*) m, itrace
! Set the break-points
dx = 2.0_nag_wp/real(nbkpts-1,kind=nag_wp)
xbkpts(1) = -1.0_nag_wp
DO i = 2, nbkpts - 1
xbkpts(i) = xbkpts(i-1) + dx
END DO
xbkpts(nbkpts) = 1.0_nag_wp
ind = 0
itask = 1
READ (nin,*) ts, tout
! Loop over output values of t
DO it = 1, 5
tout = 10.0_nag_wp*tout
! ifail: behaviour on error exit
! =0 for hard exit, =1 for quiet-soft, =-1 for noisy-soft
ifail = 0
CALL d03pdf(npde,m,ts,tout,pdedef,bndary,u,nbkpts,xbkpts,npoly,npts, &
x,uinit,acc,rsave,lrsave,isave,lisave,itask,itrace,ind,ifail)
IF (it==1) THEN
WRITE (nout,99999) npoly, nel
WRITE (nout,99998) acc, npts
WRITE (nout,99997) xout(1:6)
END IF
! Interpolate at required spatial points
ifail = 0
CALL d03pyf(npde,u,nbkpts,xbkpts,npoly,npts,xout,intpts,itype,uout, &
rsave,lrsave,ifail)
WRITE (nout,99996) ts, uout(1,1:intpts,1)
WRITE (nout,99995) uout(2,1:intpts,1)
END DO
! Print integration statistics
WRITE (nout,99994) isave(1), isave(2), isave(3), isave(5)
99999 FORMAT (' Polynomial degree =',I4,' No. of elements = ',I4)
99998 FORMAT (' Accuracy requirement =',E10.3,' Number of points = ',I5/)
99997 FORMAT (' T / X ',6F8.4/)
99996 FORMAT (1X,F7.4,' U(1)',6F8.4)
99995 FORMAT (9X,'U(2)',6F8.4/)
99994 FORMAT (' Number of integration steps in time ', &
I4/' Number of residual evaluations of resulting ODE system', &
I4/' Number of Jacobian evaluations ', &
I4/' Number of iterations of nonlinear solver ',I4)
END PROGRAM d03pdfe