fortrannumerical-methodsderivative
Large error at boudaries while calculating first and second derivative using Finite Difference
I am calculating first and second derivatives of a function a(x,y,z) stored in a 3d array a(n1,n2,n3) using Finite Difference. Here the functional value at boundaries is zero. here is the code in Fortran:
implicit none
integer i1,i2,i3
integer, parameter :: n1 = 33
integer, parameter :: n2 = 33
integer, parameter :: n3 = 32
real*8 pi, a(n1,n2,n3), a2(n1,n2,n3), z(n3),x(n1),y(n2),h,a1(n1,n2,n3)
real*8 num(n1,n2,n3),deno(n1,n2,n3),diff(n1,n2,1),A0,dx,dy
pi=3.14159265358979323846d0
dx=2.0d0*pi/(n1-1)
dy=4.0d0*pi/(n2-1)
do i1=1,n1
x(i1)=-pi+(i1-1)*dx
do i2=1,n2
y(i2)=-2.0d0*pi+(i2-1)*dy
do i3=1,n3
z(i3)=(i3-1)*2.0d0*pi/n3
a(i1,i2,1)= dcos(x(i1)/2.0d0) * dcos(y(i2)/4.0d0) !input array
a1(i1,i2,1)= - 0.25d0*dcos(x(i1)/2.0d0) * dsin(y(i2)/4.0d0) !analytical expression of first order y-derivative
enddo
enddo
enddo
do i1=1,n1
do i2=1,n2
write(20,*)x(i1),y(i2),a(i1,i2,1)
enddo
enddo
call d1y(n1,n2,n3,a,a2)
do i1=1,n1
do i2=1,n2
num(i1,i2,1)=(a2(i1,i2,1)-a1(i1,i2,1)) !numerator of error calculation
deno(i1,i2,1)=a2(i1,i2,1) !denomenator of error calculation
if (dabs(deno(i1,i2,1)) .lt. 1e-10)deno(i1,i2,1)=1.0d0
diff(i1,i2,1)=dabs(num(i1,i2,1))/dabs(deno(i1,i2,1)) !relative error in 1st order derivative calculation
write(21,*)x(i1),y(i2),a(i1,i2,1),a2(i1,i2,1),diff(i1,i2,1),a1(i1,i2,1)
write(21,*)
enddo
enddo
end
subroutine d1y(n1,n2,n3,a,a2)
implicit none
integer n1, n2, n3, i1, i2, i3
real*8 pi, a(n1,n2,n3), a2(n1,n2,n3), z(n3),x(n1),y(n2),h,a1(n1,n2,n3)
pi=3.14159265358979323846d0
h=4.0d0*pi/(n2-1)
do i1=1,n1
do i3=1,n3
do i2=1,n2
if(i2 .eq. 1)then
a2(i1,i2,i3)=( -3.0d0*a(i1,i2,i3) + 4.0d0*a(i1,i2+1,i3) - a(i1,i2+2,i3) )/ (2.0d0*h)
else if(i2 .eq. n2)then
a2(i1,i2,i3)=( 3.0d0*a(i1,i2,i3) - 4.0d0*a(i1,i2-1,i3) + a(i1,i2-2,i3) )/ (2.0d0*h)
else
a2(i1,i2,i3)=( a(i1,i2+1,i3) - a(i1,i2-1,i3) )/ (2.0d0*h)
endif
enddo
enddo
enddo
end subroutine
My input function a(i1,i2,1)= dcos(x(i1)/2.0d0) * dcos(y(i2)/4.0d0), so the sample input data (for grid no 17*17*16)
-3.1415926535897931 -6.2831853071795862 3.7493994566546440E-033
-3.1415926535897931 -5.4977871437821380 1.1945836920083898E-017
-3.1415926535897931 -4.7123889803846897 2.3432602026631496E-017
-3.1415926535897931 -3.9269908169872414 3.4018865378450254E-017
-3.1415926535897931 -3.1415926535897931 4.3297802811774670E-017
-3.1415926535897931 -2.3561944901923448 5.0912829964730140E-017
-3.1415926535897931 -1.5707963267948966 5.6571305614385013E-017
-3.1415926535897931 -0.78539816339744828 6.0055777714832775E-017
-3.1415926535897931 0.0000000000000000 6.1232339957367660E-017
-3.1415926535897931 0.78539816339744828 6.0055777714832775E-017
-3.1415926535897931 1.5707963267948966 5.6571305614385013E-017
-3.1415926535897931 2.3561944901923439 5.0912829964730146E-017
-3.1415926535897931 3.1415926535897931 4.3297802811774670E-017
-3.1415926535897931 3.9269908169872423 3.4018865378450242E-017
-3.1415926535897931 4.7123889803846897 2.3432602026631496E-017
-3.1415926535897931 5.4977871437821371 1.1945836920083910E-017
-3.1415926535897931 6.2831853071795862 3.7493994566546440E-033
-2.7488935718910690 -6.2831853071795862 1.1945836920083898E-017
-2.7488935718910690 -5.4977871437821380 3.8060233744356645E-002
-2.7488935718910690 -4.7123889803846897 7.4657834050342639E-002
-2.7488935718910690 -3.9269908169872414 0.10838637566236967
-2.7488935718910690 -3.1415926535897931 0.13794968964147156
-2.7488935718910690 -2.3561944901923448 0.16221167441072892
-2.7488935718910690 -1.5707963267948966 0.18023995550173702
-2.7488935718910690 -0.78539816339744828 0.19134171618254495
-2.7488935718910690 0.0000000000000000 0.19509032201612833
-2.7488935718910690 0.78539816339744828 0.19134171618254495
-2.7488935718910690 1.5707963267948966 0.18023995550173702
-2.7488935718910690 2.3561944901923439 0.16221167441072895
-2.7488935718910690 3.1415926535897931 0.13794968964147156
-2.7488935718910690 3.9269908169872423 0.10838637566236962
-2.7488935718910690 4.7123889803846897 7.4657834050342639E-002
-2.7488935718910690 5.4977871437821371 3.8060233744356686E-002
-2.7488935718910690 6.2831853071795862 1.1945836920083898E-017
-2.3561944901923448 -6.2831853071795862 2.3432602026631496E-017
-2.3561944901923448 -5.4977871437821380 7.4657834050342639E-002
-2.3561944901923448 -4.7123889803846897 0.14644660940672630
-2.3561944901923448 -3.9269908169872414 0.21260752369181418
-2.3561944901923448 -3.1415926535897931 0.27059805007309856
-2.3561944901923448 -2.3561944901923448 0.31818964514320852
-2.3561944901923448 -1.5707963267948966 0.35355339059327384
-2.3561944901923448 -0.78539816339744828 0.37533027751786530
-2.3561944901923448 0.0000000000000000 0.38268343236508984
-2.3561944901923448 0.78539816339744828 0.37533027751786530
-2.3561944901923448 1.5707963267948966 0.35355339059327384
-2.3561944901923448 2.3561944901923439 0.31818964514320858
-2.3561944901923448 3.1415926535897931 0.27059805007309856
-2.3561944901923448 3.9269908169872423 0.21260752369181410
-2.3561944901923448 4.7123889803846897 0.14644660940672630
-2.3561944901923448 5.4977871437821371 7.4657834050342722E-002
-2.3561944901923448 6.2831853071795862 2.3432602026631496E-017
-1.9634954084936207 -6.2831853071795862 3.4018865378450254E-017
-1.9634954084936207 -5.4977871437821380 0.10838637566236967
-1.9634954084936207 -4.7123889803846897 0.21260752369181418
-1.9634954084936207 -3.9269908169872414 0.30865828381745519
-1.9634954084936207 -3.1415926535897931 0.39284747919355117
-1.9634954084936207 -2.3561944901923448 0.46193976625564342
-1.9634954084936207 -1.5707963267948966 0.51327996715933677
-1.9634954084936207 -0.78539816339744828 0.54489510677581865
-1.9634954084936207 0.0000000000000000 0.55557023301960229
-1.9634954084936207 0.78539816339744828 0.54489510677581865
-1.9634954084936207 1.5707963267948966 0.51327996715933677
-1.9634954084936207 2.3561944901923439 0.46193976625564348
-1.9634954084936207 3.1415926535897931 0.39284747919355117
-1.9634954084936207 3.9269908169872423 0.30865828381745508
-1.9634954084936207 4.7123889803846897 0.21260752369181418
-1.9634954084936207 5.4977871437821371 0.10838637566236978
-1.9634954084936207 6.2831853071795862 3.4018865378450254E-017
-1.5707963267948966 -6.2831853071795862 4.3297802811774670E-017
-1.5707963267948966 -5.4977871437821380 0.13794968964147156
-1.5707963267948966 -4.7123889803846897 0.27059805007309856
-1.5707963267948966 -3.9269908169872414 0.39284747919355117
-1.5707963267948966 -3.1415926535897931 0.50000000000000011
-1.5707963267948966 -2.3561944901923448 0.58793780120967942
-1.5707963267948966 -1.5707963267948966 0.65328148243818829
-1.5707963267948966 -0.78539816339744828 0.69351992266107376
-1.5707963267948966 0.0000000000000000 0.70710678118654757
-1.5707963267948966 0.78539816339744828 0.69351992266107376
-1.5707963267948966 1.5707963267948966 0.65328148243818829
-1.5707963267948966 2.3561944901923439 0.58793780120967942
-1.5707963267948966 3.1415926535897931 0.50000000000000011
-1.5707963267948966 3.9269908169872423 0.39284747919355101
-1.5707963267948966 4.7123889803846897 0.27059805007309856
-1.5707963267948966 5.4977871437821371 0.13794968964147170
-1.5707963267948966 6.2831853071795862 4.3297802811774670E-017
-1.1780972450961724 -6.2831853071795862 5.0912829964730140E-017
-1.1780972450961724 -5.4977871437821380 0.16221167441072892
-1.1780972450961724 -4.7123889803846897 0.31818964514320852
-1.1780972450961724 -3.9269908169872414 0.46193976625564342
-1.1780972450961724 -3.1415926535897931 0.58793780120967942
-1.1780972450961724 -2.3561944901923448 0.69134171618254492
-1.1780972450961724 -1.5707963267948966 0.76817775671141630
-1.1780972450961724 -0.78539816339744828 0.81549315684891710
-1.1780972450961724 0.0000000000000000 0.83146961230254524
-1.1780972450961724 0.78539816339744828 0.81549315684891710
-1.1780972450961724 1.5707963267948966 0.76817775671141630
-1.1780972450961724 2.3561944901923439 0.69134171618254503
-1.1780972450961724 3.1415926535897931 0.58793780120967942
-1.1780972450961724 3.9269908169872423 0.46193976625564326
-1.1780972450961724 4.7123889803846897 0.31818964514320852
-1.1780972450961724 5.4977871437821371 0.16221167441072909
-1.1780972450961724 6.2831853071795862 5.0912829964730140E-017
-0.78539816339744828 -6.2831853071795862 5.6571305614385013E-017
-0.78539816339744828 -5.4977871437821380 0.18023995550173702
-0.78539816339744828 -4.7123889803846897 0.35355339059327384
-0.78539816339744828 -3.9269908169872414 0.51327996715933677
-0.78539816339744828 -3.1415926535897931 0.65328148243818829
-0.78539816339744828 -2.3561944901923448 0.76817775671141630
-0.78539816339744828 -1.5707963267948966 0.85355339059327373
-0.78539816339744828 -0.78539816339744828 0.90612744635288778
-0.78539816339744828 0.0000000000000000 0.92387953251128674
-0.78539816339744828 0.78539816339744828 0.90612744635288778
-0.78539816339744828 1.5707963267948966 0.85355339059327373
-0.78539816339744828 2.3561944901923439 0.76817775671141642
-0.78539816339744828 3.1415926535897931 0.65328148243818829
-0.78539816339744828 3.9269908169872423 0.51327996715933655
-0.78539816339744828 4.7123889803846897 0.35355339059327384
-0.78539816339744828 5.4977871437821371 0.18023995550173721
-0.78539816339744828 6.2831853071795862 5.6571305614385013E-017
-0.39269908169872414 -6.2831853071795862 6.0055777714832775E-017
-0.39269908169872414 -5.4977871437821380 0.19134171618254495
-0.39269908169872414 -4.7123889803846897 0.37533027751786530
-0.39269908169872414 -3.9269908169872414 0.54489510677581865
-0.39269908169872414 -3.1415926535897931 0.69351992266107376
-0.39269908169872414 -2.3561944901923448 0.81549315684891710
-0.39269908169872414 -1.5707963267948966 0.90612744635288778
-0.39269908169872414 -0.78539816339744828 0.96193976625564337
-0.39269908169872414 0.0000000000000000 0.98078528040323043
-0.39269908169872414 0.78539816339744828 0.96193976625564337
-0.39269908169872414 1.5707963267948966 0.90612744635288778
-0.39269908169872414 2.3561944901923439 0.81549315684891721
-0.39269908169872414 3.1415926535897931 0.69351992266107376
-0.39269908169872414 3.9269908169872423 0.54489510677581843
-0.39269908169872414 4.7123889803846897 0.37533027751786530
-0.39269908169872414 5.4977871437821371 0.19134171618254514
-0.39269908169872414 6.2831853071795862 6.0055777714832775E-017
0.0000000000000000 -6.2831853071795862 6.1232339957367660E-017
0.0000000000000000 -5.4977871437821380 0.19509032201612833
0.0000000000000000 -4.7123889803846897 0.38268343236508984
0.0000000000000000 -3.9269908169872414 0.55557023301960229
0.0000000000000000 -3.1415926535897931 0.70710678118654757
0.0000000000000000 -2.3561944901923448 0.83146961230254524
0.0000000000000000 -1.5707963267948966 0.92387953251128674
0.0000000000000000 -0.78539816339744828 0.98078528040323043
0.0000000000000000 0.0000000000000000 1.0000000000000000
0.0000000000000000 0.78539816339744828 0.98078528040323043
0.0000000000000000 1.5707963267948966 0.92387953251128674
0.0000000000000000 2.3561944901923439 0.83146961230254535
0.0000000000000000 3.1415926535897931 0.70710678118654757
0.0000000000000000 3.9269908169872423 0.55557023301960207
0.0000000000000000 4.7123889803846897 0.38268343236508984
0.0000000000000000 5.4977871437821371 0.19509032201612853
0.0000000000000000 6.2831853071795862 6.1232339957367660E-017
0.39269908169872414 -6.2831853071795862 6.0055777714832775E-017
0.39269908169872414 -5.4977871437821380 0.19134171618254495
0.39269908169872414 -4.7123889803846897 0.37533027751786530
0.39269908169872414 -3.9269908169872414 0.54489510677581865
0.39269908169872414 -3.1415926535897931 0.69351992266107376
0.39269908169872414 -2.3561944901923448 0.81549315684891710
0.39269908169872414 -1.5707963267948966 0.90612744635288778
0.39269908169872414 -0.78539816339744828 0.96193976625564337
0.39269908169872414 0.0000000000000000 0.98078528040323043
0.39269908169872414 0.78539816339744828 0.96193976625564337
0.39269908169872414 1.5707963267948966 0.90612744635288778
0.39269908169872414 2.3561944901923439 0.81549315684891721
0.39269908169872414 3.1415926535897931 0.69351992266107376
0.39269908169872414 3.9269908169872423 0.54489510677581843
0.39269908169872414 4.7123889803846897 0.37533027751786530
0.39269908169872414 5.4977871437821371 0.19134171618254514
0.39269908169872414 6.2831853071795862 6.0055777714832775E-017
0.78539816339744828 -6.2831853071795862 5.6571305614385013E-017
0.78539816339744828 -5.4977871437821380 0.18023995550173702
0.78539816339744828 -4.7123889803846897 0.35355339059327384
0.78539816339744828 -3.9269908169872414 0.51327996715933677
0.78539816339744828 -3.1415926535897931 0.65328148243818829
0.78539816339744828 -2.3561944901923448 0.76817775671141630
0.78539816339744828 -1.5707963267948966 0.85355339059327373
0.78539816339744828 -0.78539816339744828 0.90612744635288778
0.78539816339744828 0.0000000000000000 0.92387953251128674
0.78539816339744828 0.78539816339744828 0.90612744635288778
0.78539816339744828 1.5707963267948966 0.85355339059327373
0.78539816339744828 2.3561944901923439 0.76817775671141642
0.78539816339744828 3.1415926535897931 0.65328148243818829
0.78539816339744828 3.9269908169872423 0.51327996715933655
0.78539816339744828 4.7123889803846897 0.35355339059327384
0.78539816339744828 5.4977871437821371 0.18023995550173721
0.78539816339744828 6.2831853071795862 5.6571305614385013E-017
1.1780972450961720 -6.2831853071795862 5.0912829964730146E-017
1.1780972450961720 -5.4977871437821380 0.16221167441072895
1.1780972450961720 -4.7123889803846897 0.31818964514320858
1.1780972450961720 -3.9269908169872414 0.46193976625564348
1.1780972450961720 -3.1415926535897931 0.58793780120967942
1.1780972450961720 -2.3561944901923448 0.69134171618254503
1.1780972450961720 -1.5707963267948966 0.76817775671141642
1.1780972450961720 -0.78539816339744828 0.81549315684891721
1.1780972450961720 0.0000000000000000 0.83146961230254535
1.1780972450961720 0.78539816339744828 0.81549315684891721
1.1780972450961720 1.5707963267948966 0.76817775671141642
1.1780972450961720 2.3561944901923439 0.69134171618254503
1.1780972450961720 3.1415926535897931 0.58793780120967942
1.1780972450961720 3.9269908169872423 0.46193976625564331
1.1780972450961720 4.7123889803846897 0.31818964514320858
1.1780972450961720 5.4977871437821371 0.16221167441072912
1.1780972450961720 6.2831853071795862 5.0912829964730146E-017
1.5707963267948966 -6.2831853071795862 4.3297802811774670E-017
1.5707963267948966 -5.4977871437821380 0.13794968964147156
1.5707963267948966 -4.7123889803846897 0.27059805007309856
1.5707963267948966 -3.9269908169872414 0.39284747919355117
1.5707963267948966 -3.1415926535897931 0.50000000000000011
1.5707963267948966 -2.3561944901923448 0.58793780120967942
1.5707963267948966 -1.5707963267948966 0.65328148243818829
1.5707963267948966 -0.78539816339744828 0.69351992266107376
1.5707963267948966 0.0000000000000000 0.70710678118654757
1.5707963267948966 0.78539816339744828 0.69351992266107376
1.5707963267948966 1.5707963267948966 0.65328148243818829
1.5707963267948966 2.3561944901923439 0.58793780120967942
1.5707963267948966 3.1415926535897931 0.50000000000000011
1.5707963267948966 3.9269908169872423 0.39284747919355101
1.5707963267948966 4.7123889803846897 0.27059805007309856
1.5707963267948966 5.4977871437821371 0.13794968964147170
1.5707963267948966 6.2831853071795862 4.3297802811774670E-017
1.9634954084936211 -6.2831853071795862 3.4018865378450242E-017
1.9634954084936211 -5.4977871437821380 0.10838637566236962
1.9634954084936211 -4.7123889803846897 0.21260752369181410
1.9634954084936211 -3.9269908169872414 0.30865828381745508
1.9634954084936211 -3.1415926535897931 0.39284747919355101
1.9634954084936211 -2.3561944901923448 0.46193976625564326
1.9634954084936211 -1.5707963267948966 0.51327996715933655
1.9634954084936211 -0.78539816339744828 0.54489510677581843
1.9634954084936211 0.0000000000000000 0.55557023301960207
1.9634954084936211 0.78539816339744828 0.54489510677581843
1.9634954084936211 1.5707963267948966 0.51327996715933655
1.9634954084936211 2.3561944901923439 0.46193976625564331
1.9634954084936211 3.1415926535897931 0.39284747919355101
1.9634954084936211 3.9269908169872423 0.30865828381745491
1.9634954084936211 4.7123889803846897 0.21260752369181410
1.9634954084936211 5.4977871437821371 0.10838637566236972
1.9634954084936211 6.2831853071795862 3.4018865378450242E-017
2.3561944901923448 -6.2831853071795862 2.3432602026631496E-017
2.3561944901923448 -5.4977871437821380 7.4657834050342639E-002
2.3561944901923448 -4.7123889803846897 0.14644660940672630
2.3561944901923448 -3.9269908169872414 0.21260752369181418
2.3561944901923448 -3.1415926535897931 0.27059805007309856
2.3561944901923448 -2.3561944901923448 0.31818964514320852
2.3561944901923448 -1.5707963267948966 0.35355339059327384
2.3561944901923448 -0.78539816339744828 0.37533027751786530
2.3561944901923448 0.0000000000000000 0.38268343236508984
2.3561944901923448 0.78539816339744828 0.37533027751786530
2.3561944901923448 1.5707963267948966 0.35355339059327384
2.3561944901923448 2.3561944901923439 0.31818964514320858
2.3561944901923448 3.1415926535897931 0.27059805007309856
2.3561944901923448 3.9269908169872423 0.21260752369181410
2.3561944901923448 4.7123889803846897 0.14644660940672630
2.3561944901923448 5.4977871437821371 7.4657834050342722E-002
2.3561944901923448 6.2831853071795862 2.3432602026631496E-017
2.7488935718910685 -6.2831853071795862 1.1945836920083910E-017
2.7488935718910685 -5.4977871437821380 3.8060233744356686E-002
2.7488935718910685 -4.7123889803846897 7.4657834050342722E-002
2.7488935718910685 -3.9269908169872414 0.10838637566236978
2.7488935718910685 -3.1415926535897931 0.13794968964147170
2.7488935718910685 -2.3561944901923448 0.16221167441072909
2.7488935718910685 -1.5707963267948966 0.18023995550173721
2.7488935718910685 -0.78539816339744828 0.19134171618254514
2.7488935718910685 0.0000000000000000 0.19509032201612853
2.7488935718910685 0.78539816339744828 0.19134171618254514
2.7488935718910685 1.5707963267948966 0.18023995550173721
2.7488935718910685 2.3561944901923439 0.16221167441072912
2.7488935718910685 3.1415926535897931 0.13794968964147170
2.7488935718910685 3.9269908169872423 0.10838637566236972
2.7488935718910685 4.7123889803846897 7.4657834050342722E-002
2.7488935718910685 5.4977871437821371 3.8060233744356721E-002
2.7488935718910685 6.2831853071795862 1.1945836920083910E-017
3.1415926535897931 -6.2831853071795862 3.7493994566546440E-033
3.1415926535897931 -5.4977871437821380 1.1945836920083898E-017
3.1415926535897931 -4.7123889803846897 2.3432602026631496E-017
3.1415926535897931 -3.9269908169872414 3.4018865378450254E-017
3.1415926535897931 -3.1415926535897931 4.3297802811774670E-017
3.1415926535897931 -2.3561944901923448 5.0912829964730140E-017
3.1415926535897931 -1.5707963267948966 5.6571305614385013E-017
3.1415926535897931 -0.78539816339744828 6.0055777714832775E-017
3.1415926535897931 0.0000000000000000 6.1232339957367660E-017
3.1415926535897931 0.78539816339744828 6.0055777714832775E-017
3.1415926535897931 1.5707963267948966 5.6571305614385013E-017
3.1415926535897931 2.3561944901923439 5.0912829964730146E-017
3.1415926535897931 3.1415926535897931 4.3297802811774670E-017
3.1415926535897931 3.9269908169872423 3.4018865378450242E-017
3.1415926535897931 4.7123889803846897 2.3432602026631496E-017
3.1415926535897931 5.4977871437821371 1.1945836920083910E-017
3.1415926535897931 6.2831853071795862 3.7493994566546440E-033
Input and output functions are as shown in the figure below.
As the output function is the same as - 0.25d0*dcos(x(i1)/2.0d0) * dsin(y(i2)/4.0d0), this shows derivative calculation is correct. Here in subroutine 'd1y', I have used forward and backward finite difference formula for calculation of derivative at boundaries and the central difference for points lying between two boundaries. Then I have calculated the relative error. I got an error of 0.003 at the boundary of y axis and 0.0015 at the points in between boundaries as shown in fig below
.
I calculated 2nd derivative using same technique. The subroutine is given below:
`subroutine d2y(n1,n2,n3,a,a2)
implicit none
integer n1, n2, n3, i1, i2, i3
real*8 pi, a(n1,n2,n3), a2(n1,n2,n3), z(n3),x(n1),y(n2),h
pi=3.14159265358979323846d0
h=4.0d0*pi/(n2-1)
a2=0.0d0
i3 =1
do i1=1,n1
do i2=1,n2
if(i2 == 1)then
a2(i1,i2,i3)=( 2.0d0*a(i1,i2,i3) - 5.0d0*a(i1,i2+1,i3) + 4.0d0*a(i1,i2+2,i3) - a(i1,i2+3,i3))/(h*h)
else if( i2 == n2)then
a2(i1,i2,i3)=( 2.0d0*a(i1,i2,i3) - 5.0d0*a(i1,i2-1,i3) + 4.0d0*a(i1,i2-2,i3) - a(i1,i2-3,i3))/(h*h)
else
a2(i1,i2,i3)= ( a(i1,i2+1,i3) - 2.0d0*a(i1,i2,i3) + a(i1,i2-1,i3) )/(h*h)
endif
enddo
enddo
! enddo
end subroutine
Here also the error is large at boundaries, in fact, the error is very large compared to first order derivative as shown in figure:
.
Why this large error in boundary? Please explain.
`
Solution
First, I stand corrected. The names 'Forward Difference' and 'Backward Difference' usually refer to the standard 1st order formulas. You do have correct 2nd order, One-sided 1st derivative formulas.
The 2nd order, 1st derivative, Forward difference formula is derived with a Taylor series expression at x+h and at x+2h as follows:
f(x+h) = f(x) + h*f'(x) + h^2*f''(x)/2! + h^3*f'''(x)/3! ....
f(x+2h) = f(x) + 2h*f'(x) + 4h^2*f''(x)/2! + 8h^3*f'''(x)/3! ....
Now, take 4 times the 1st series and subtract the 2nd series and solve for f'(x). Note that the second derivative term disappears.
f'(x) = 3h/2 * f(x) - 2/h *f(x+h) + 1/(2h) * f(x+2h) + 0 - 4h^2*f'''(x)/3!
This is the formula you are using and it is 2nd order accurate. This does not mean that the accuracy is identical to the central difference accuracy, only that the order of the accuracy is the same.
If you look at the derivation of central difference while keeping the error term, you will find that the error term looks like:
-h^2*f'''(x)/6
Did you note that the error term in the forward difference formula has a constant in front of it. In general, the error will always be larger. But that is not what it means to be 2nd order accurate. The order of accuracy tells you how the error changes when you change h.
For example, if you halve h, you should get about 4 times the accuracy. That means in your example above, the central difference error will go from about 0.0015 to about 0.0004 and the boundary error will go from about 0.003 to about 0.0008
Final Takeaway: Your program is correct and your errors are correct!
- Can two FORTRAN modules use parts of one another?
- Changing derived type data members in `SUBROUTINE` is allowed even with INTENT(IN). How to protect data?
- malloc(): unaligned tcache chunk detected. Has anyone faced this before for MPI fortran programs?
- Function Returning an array in Fortran
- Is there a limit to the number of arguments passed to a fortran function?
- Pass a fortran derived type which contains allocatable array between different compilers(PGI and Intel)
- Numerical discrepancies after upgrading Fortran compilers (Intel v11.1 to v13.0)
- Correct suffix for Fortran 2003 source file - Intel Fortran compiler
- Cannot execute wfc386 no such file or directory - Watcom Fortran IDE
- Should casting from ``real`` to ``double precision`` in FORTRAN be affected by the compiler?
- "ftn" command in ubuntu system
- How can I force an array overflow?
- Computing the cross product of two vectors in Fortran 90
- PGI Flang github disappeared
- gfortran Deallocation segmentation fault
- How to avoid skipping input with Fortran READ statement
- Best way to handle multiple exit points in Fortran
- impure function inside a DO CONCURRENT block
- Intel Compilers: linking a fortran compiled library to main program of C
- How to automatically adjust string length in the A edit descriptor?
- Is it safe to pass a pointer of a non-contiguous array to a Lapack95 routine in Fortran?
- How to print a colored text without including ANSI escape sequence in the output length?
- Is it possible to somehow mix static and dynamic loop scheduling?
- How to compile Fortran 77 files in linux?
- Compiler validation suite for Fortran 2003?
- gfortran performs unsafe math optimizations of cos(atan(x))
- Polymorphic assumed-rank entity used in SELECT TYPE
- What are TRUE and FALSE constants (without the surrounding periods ('.')) in Fortran?
- Compile a library that uses mpi (AGMG)
- What does 10** mean in Fortran?