Computational Physics: Fortran Version, 计算物理的经典之作。认认真真的读这本书,好好的做笔记。读本书第一个example(长案例)便让人获益匪浅。让我瞬间感觉到科班和野路子就是不一样啊!
书中的代码是77版的,读起来有点不爽。试着改写成95版的。
Chapter 1 Basic Mathematical Operations
求导数,求面积(与积分这个词略有差异,至少本书里是这么约定的),寻根:三种很基本的要素。
这三种要素是大部分计算模型的核心。
1.1 求导数(差分)
(1)3点公式
(2)误差更大些的“2点公式”:前向差分公式,后向差分公式
(3)5点公式
下面是几个公式的比较,书中有77版代码。
-
!ccccccccccccccccccccccccccccccccccccccccccccccccc!
-
! chap1a.f90
-
! Computational Physics: Fortran Version
-
! 3-point formula,2-point formulas,5-point formula
-
program main
-
!ccccccccccccccccccccccccccccccccccccccccccccccccc!
-
!
-
!variables:
-
implicit none
-
integer :: counter ! using in do loop
-
integer,parameter :: num=100 ! using in do loop
-
real :: x=1.0,exact ! x: variable. exact: cos(1.0).
-
real :: threeP,twoFP,twoBP,fiveP ! threeP: three-point. twoFP: two-point forward.
-
! twoBP: two-point backward. fiveP: five-point
-
real :: diff1,diff2,diff3,diff4 ! diff=cos(1.0)-approximate
-
real :: h=0.1 ! the step in calculation
-
!
-
exact=cos(x)
-
!
-
do counter=1,num,1
-
threeP=(sin(x+h)-sin(x-h))/(2*h)
-
diff1=exact-threeP
-
write(*,*)"exact:",exact,"threeP:",threeP,"diff1:",diff1
-
!
-
twoFP=(sin(x+h)-sin(x))/h
-
diff2=exact-twoFP
-
write(*,*)"exact:",exact,"twoFP:",twoFP,"diff2:",diff2
-
!
-
twoBP=(sin(x)-sin(x-h))/h
-
diff3=exact-twoBP
-
write(*,*)"exact:",exact,"twoBP:",twoBP,"diff3:",diff3
-
!
-
fiveP=(sin(x-2*h)-8*sin(x-h)+8*sin(x+h)-sin(x+2*h))/(12*h)
-
diff4=exact-fiveP
-
write(*,*)"exact:",exact,"fiveP:",fiveP,"diff4:",diff4
-
h=h/2.0
-
write(*,*)"*****************","over of the",counter,"th line","**********************************"
-
!
-
end do
-
!
-
stop
-
-
end program main
输出结果的前20行,误差已经经历过一个戏剧性的变化:由大变小,再次变得更大,最后大到无法忍受。
exact: 0.540302277 threeP: 0.539402485 diff1: 8.99791718E-04
exact: 0.540302277 twoFP: 0.497364104 diff2: 4.29381728E-02
exact: 0.540302277 twoBP: 0.581440628 diff3: -4.11383510E-02
exact: 0.540302277 fiveP: 0.540300548 diff4: 1.72853470E-06
*****************over of the 1 th line**********************************
exact: 0.540302277 threeP: 0.540076792 diff1: 2.25484371E-04
exact: 0.540302277 twoFP: 0.519043803 diff2: 2.12584734E-02
exact: 0.540302277 twoBP: 0.561109185 diff3: -2.08069086E-02
exact: 0.540302277 fiveP: 0.540301919 diff4: 3.57627869E-07
*****************over of the 2 th line**********************************
exact: 0.540302277 threeP: 0.540245056 diff1: 5.72204590E-05
exact: 0.540302277 twoFP: 0.529726803 diff2: 1.05754733E-02
exact: 0.540302277 twoBP: 0.550762177 diff3: -1.04598999E-02
exact: 0.540302277 fiveP: 0.540300906 diff4: 1.37090683E-06
*****************over of the 3 th line**********************************
exact: 0.540302277 threeP: 0.540289879 diff1: 1.23977661E-05
exact: 0.540302277 twoFP: 0.535032868 diff2: 5.26940823E-03
exact: 0.540302277 twoBP: 0.545544565 diff3: -5.24228811E-03
exact: 0.540302277 fiveP: 0.540304422 diff4: -2.14576721E-06
*****************over of the 4 th line**********************************
exact: 0.540302277 threeP: 0.540298998 diff1: 3.27825546E-06
exact: 0.540302277 twoFP: 0.537667572 diff2: 2.63470411E-03
exact: 0.540302277 twoBP: 0.542925954 diff3: -2.62367725E-03
exact: 0.540302277 fiveP: 0.540302873 diff4: -5.96046448E-07
*****************over of the 5 th line**********************************
exact: 0.540302277 threeP: 0.540293932 diff1: 8.34465027E-06
exact: 0.540302277 twoFP: 0.538969636 diff2: 1.33264065E-03
exact: 0.540302277 twoBP: 0.541609347 diff3: -1.30707026E-03
exact: 0.540302277 fiveP: 0.540295601 diff4: 6.67572021E-06
*****************over of the 6 th line**********************************
exact: 0.540302277 threeP: 0.540299237 diff1: 3.03983688E-06
exact: 0.540302277 twoFP: 0.539647281 diff2: 6.54995441E-04
exact: 0.540302277 twoBP: 0.540933311 diff3: -6.31034374E-04
exact: 0.540302277 fiveP: 0.540286839 diff4: 1.54376030E-05
*****************over of the 7 th line**********************************
exact: 0.540302277 threeP: 0.540297687 diff1: 4.58955765E-06
exact: 0.540302277 twoFP: 0.539972663 diff2: 3.29613686E-04
exact: 0.540302277 twoBP: 0.540586829 diff3: -2.84552574E-04
exact: 0.540302277 fiveP: 0.540347278 diff4: -4.50015068E-05
*****************over of the 8 th line**********************************
exact: 0.540302277 threeP: 0.540370822 diff1: -6.85453415E-05
exact: 0.540302277 twoFP: 0.540242016 diff2: 6.02602959E-05
exact: 0.540302277 twoBP: 0.540427923 diff3: -1.25646591E-04
exact: 0.540302277 fiveP: 0.540282667 diff4: 1.96099281E-05
*****************over of the 9 th line**********************************
exact: 0.540302277 threeP: 0.540217578 diff1: 8.46982002E-05
exact: 0.540302277 twoFP: 0.540017724 diff2: 2.84552574E-04
exact: 0.540302277 twoBP: 0.540274084 diff3: 2.81929970E-05
exact: 0.540302277 fiveP: 0.540198386 diff4: 1.03890896E-04
*****************over of the 10 th line**********************************
exact: 0.540302277 threeP: 0.540042877 diff1: 2.59399414E-04
exact: 0.540302277 twoFP: 0.539874375 diff2: 4.27901745E-04
exact: 0.540302277 twoBP: 0.539924681 diff3: 3.77595425E-04
exact: 0.540302277 fiveP: 0.540200233 diff4: 1.02043152E-04
*****************over of the 11 th line**********************************
exact: 0.540302277 threeP: 0.540804744 diff1: -5.02467155E-04
exact: 0.540302277 twoFP: 0.541418612 diff2: -1.11633539E-03
exact: 0.540302277 twoBP: 0.539617360 diff3: 6.84916973E-04
exact: 0.540302277 fiveP: 0.540627241 diff4: -3.24964523E-04
*****************over of the 12 th line**********************************
exact: 0.540302277 threeP: 0.540842652 diff1: -5.40375710E-04
exact: 0.540302277 twoFP: 0.540845096 diff2: -5.42819500E-04
exact: 0.540302277 twoBP: 0.539693117 diff3: 6.09159470E-04
exact: 0.540302277 fiveP: 0.540466130 diff4: -1.63853168E-04
*****************over of the 13 th line**********************************
exact: 0.540302277 threeP: 0.540266514 diff1: 3.57627869E-05
exact: 0.540302277 twoFP: 0.539697945 diff2: 6.04331493E-04
exact: 0.540302277 twoBP: 0.538540840 diff3: 1.76143646E-03
exact: 0.540302277 fiveP: 0.538741231 diff4: 1.56104565E-03
*****************over of the 14 th line**********************************
exact: 0.540302277 threeP: 0.540240407 diff1: 6.18696213E-05
exact: 0.540302277 twoFP: 0.542286575 diff2: -1.98429823E-03
exact: 0.540302277 twoBP: 0.533605874 diff3: 6.69640303E-03
exact: 0.540302277 fiveP: 0.535733640 diff4: 4.56863642E-03
*****************over of the 15 th line**********************************
exact: 0.540302277 threeP: 0.547711194 diff1: -7.40891695E-03
exact: 0.540302277 twoFP: 0.557229400 diff2: -1.69271231E-02
exact: 0.540302277 twoBP: 0.529016137 diff3: 1.12861395E-02
exact: 0.540302277 fiveP: 0.549185514 diff4: -8.88323784E-03
*****************over of the 16 th line**********************************
exact: 0.540302277 threeP: 0.558164477 diff1: -1.78622007E-02
exact: 0.540302277 twoFP: 0.567583740 diff2: -2.72814631E-02
exact: 0.540302277 twoBP: 0.530391455 diff3: 9.91082191E-03
exact: 0.540302277 fiveP: 0.558054984 diff4: -1.77527070E-02
*****************over of the 17 th line**********************************
exact: 0.540302277 threeP: 0.529456198 diff1: 1.08460784E-02
exact: 0.540302277 twoFP: 0.510167480 diff2: 3.01347971E-02
exact: 0.540302277 twoBP: 0.512037396 diff3: 2.82648802E-02
exact: 0.540302277 fiveP: 0.504442036 diff4: 3.58602405E-02
*****************over of the 18 th line**********************************
exact: 0.540302277 threeP: 0.529059291 diff1: 1.12429857E-02
exact: 0.540302277 twoFP: 0.551585019 diff2: -1.12827420E-02
exact: 0.540302277 twoBP: 0.433118582 diff3: 0.107183695
exact: 0.540302277 fiveP: 0.456334889 diff4: 8.39673877E-02
*****************over of the 19 th line**********************************
exact: 0.540302277 threeP: 0.648601770 diff1: -0.108299494
exact: 0.540302277 twoFP: 0.790670037 diff2: -0.250367761
exact: 0.540302277 twoBP: 0.359703541 diff3: 0.180598736
exact: 0.540302277 fiveP: 0.554383636 diff4: -1.40813589E-02
*****************over of the 20 th line**********************************
exact: 0.540302277 threeP: 0.815859020 diff1: -0.275556743
exact: 0.540302277 twoFP: 0.956340134 diff2: -0.416037858
exact: 0.540302277 twoBP: 0.381718069 diff3: 0.158584207
exact: 0.540302277 fiveP: 0.696298718 diff4: -0.155996442
改变变量的存储长度,kind=8(双精度可以使得数值部分达到15~16位有效数字)之后呢?一直到第50行,误差才达到以上的程度。
关于舍入误差以及其他误差来源,以后再专门学习一下吧。