使用有限差分计算一阶和二阶导数时边界处存在较大误差
Large error at boudaries while calculating first and second derivative using Finite Difference
我正在使用有限差分计算存储在 3d 数组 a(n1,n2,n3) 中的函数 a(x,y,z) 的一阶和二阶导数。这里边界处的函数值为零。这是 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
我的输入函数a(i1,i2,1)= dcos(x(i1)/2.0d0) * dcos(y(i2)/4.0d0)
,所以示例输入数据(对于网格号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
输入输出函数如下图所示。
由于输出函数与- 0.25d0*dcos(x(i1)/2.0d0) * dsin(y(i2)/4.0d0)
相同,说明导数计算正确。在子程序'd1y'中,我使用前向和后向有限差分公式计算边界处的导数和位于两个边界之间的点的中心差。然后我计算了相对误差。我在 y 轴的边界处得到 0.003 的错误,在边界之间的点处得到 0.0015 的错误,如下图所示。
我使用相同的技术计算了二阶导数。子程序如下:
`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
这里同样在边界处误差很大,其实和一阶导数相比误差很大如图:.
为什么边界有这么大的错误?请解释。
`
首先,我的观点是正确的。名称 'Forward Difference' 和 'Backward Difference' 通常指的是标准的一阶公式。你确实有正确的二阶单边一阶导数公式。
二阶一阶导数前向差分公式在x+h和x+2h处用泰勒级数表达式推导如下:
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! ....
现在,用第一个数列的 4 倍减去第二个数列,求解 f'(x)。请注意,二阶导数项消失了。
f'(x) = 3h/2 * f(x) - 2/h *f(x+h) + 1/(2h) * f(x+2h) + 0 - 4h^2* f'''(x)/3!
这是您使用的公式,二阶准确。这并不意味着精度与中心差分精度相同,只是精度的阶数相同。
如果你在保留误差项的情况下看中心差的推导,你会发现误差项是这样的:
-h^2*f'''(x)/6
你有没有注意到前向差分公式中的误差项前面有一个常数。一般来说,误差总是会更大。但这并不是二阶准确的意思。精度的顺序告诉你当你改变 h 时误差如何变化。
例如,如果将 h 减半,则应获得大约 4 倍的准确度。这意味着在上面的示例中,中心差分误差将从大约 0.0015 变为大约 0.0004,边界误差将从大约 0.003 变为大约 0.0008
最后要点:你的程序是正确的,你的错误也是正确的!
我正在使用有限差分计算存储在 3d 数组 a(n1,n2,n3) 中的函数 a(x,y,z) 的一阶和二阶导数。这里边界处的函数值为零。这是 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
我的输入函数a(i1,i2,1)= dcos(x(i1)/2.0d0) * dcos(y(i2)/4.0d0)
,所以示例输入数据(对于网格号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
输入输出函数如下图所示。
由于输出函数与- 0.25d0*dcos(x(i1)/2.0d0) * dsin(y(i2)/4.0d0)
相同,说明导数计算正确。在子程序'd1y'中,我使用前向和后向有限差分公式计算边界处的导数和位于两个边界之间的点的中心差。然后我计算了相对误差。我在 y 轴的边界处得到 0.003 的错误,在边界之间的点处得到 0.0015 的错误,如下图所示
我使用相同的技术计算了二阶导数。子程序如下:
`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
这里同样在边界处误差很大,其实和一阶导数相比误差很大如图:
为什么边界有这么大的错误?请解释。
`
首先,我的观点是正确的。名称 'Forward Difference' 和 'Backward Difference' 通常指的是标准的一阶公式。你确实有正确的二阶单边一阶导数公式。
二阶一阶导数前向差分公式在x+h和x+2h处用泰勒级数表达式推导如下:
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! ....
现在,用第一个数列的 4 倍减去第二个数列,求解 f'(x)。请注意,二阶导数项消失了。
f'(x) = 3h/2 * f(x) - 2/h *f(x+h) + 1/(2h) * f(x+2h) + 0 - 4h^2* f'''(x)/3!
这是您使用的公式,二阶准确。这并不意味着精度与中心差分精度相同,只是精度的阶数相同。
如果你在保留误差项的情况下看中心差的推导,你会发现误差项是这样的:
-h^2*f'''(x)/6
你有没有注意到前向差分公式中的误差项前面有一个常数。一般来说,误差总是会更大。但这并不是二阶准确的意思。精度的顺序告诉你当你改变 h 时误差如何变化。
例如,如果将 h 减半,则应获得大约 4 倍的准确度。这意味着在上面的示例中,中心差分误差将从大约 0.0015 变为大约 0.0004,边界误差将从大约 0.003 变为大约 0.0008
最后要点:你的程序是正确的,你的错误也是正确的!