Keyword: チェビシェフ, 選点法, 放物型偏微分方程式
概要
本サンプルはチェビシェフC0選点法を用いた離散化により放物型偏微分方程式を解くFortranによるサンプルプログラムです。 本サンプルは以下に示される1組の2次の放物型偏微分方程式(PDE)で記述される4次の偏微分方程式(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
- 1行目はタイトル行で読み飛ばされます。
- 2行目は補間数(intpts=6)、ブレークポイントの数(nbkpts=10)、偏微分方程式(PDE)の解の近似に使用されるチェビシェフ多項式の次数(npoly=3)、実行される補間(itype=1:補間点の解を計算)を指定しています。
- 3行目は空間補間点(xout)を指定しています。
- 4行目は局所誤差を制御するための正数(acc)を指定しています。
- 5行目は座標システム(m=0:デカルト座標)とトレース情報のレベル(itrace=0:警告メッセージのみ出力)を指定しています。
- 6行目は独立変数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
- 2行目にはチェビシェフ多項式の次数と要素の数が出力されています。
- 3行目には誤差推定を制御するための正数とメッシュ点の数が出力されています。
- 5行目にはブレークポイントが出力されています。
- 7~20行目には独立変数の値とU1(x,t)とU2(x,t)の計算結果が出力されています。
- 22行目には積分のステップの数が出力されています。
- 23行目には結果として生じる常微分方程式(ODE)の残差評価の数が出力されています。
- 24行目にはヤコビアン評価の数が出力されています。
- 25行目には非線形ソルバの反復数が出力されています。
ソースコード
(本ルーチンの詳細は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