チェビシェフC0選点法を用いた離散化による放物型偏微分方程式

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

ホーム > 製品 > nAG数値計算ライブラリ > nAG Fortranライブラリ > サンプルソースコード集 > チェビシェフC0選点法を用いた離散化による放物型偏微分方程式

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


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