Chinaunix首页 | 论坛 | 博客
  • 博客访问: 40921
  • 博文数量: 10
  • 博客积分: 0
  • 博客等级: 民兵
  • 技术积分: 95
  • 用 户 组: 普通用户
  • 注册时间: 2015-10-18 15:21
文章分类

全部博文(10)

文章存档

2015年(10)

我的朋友

分类: C/C++

2015-10-24 10:47:36

         Computational Physics: Fortran Version, 计算物理的经典之作。认认真真的读这本书,好好的做笔记。读本书第一个example(长案例)便让人获益匪浅。让我瞬间感觉到科班和野路子就是不一样啊!
         书中的代码是77版的,读起来有点不爽。试着改写成95版的。
Chapter 1       Basic Mathematical Operations
  求导数,求面积(与积分这个词略有差异,至少本书里是这么约定的),寻根:三种很基本的要素。
  这三种要素是大部分计算模型的核心。
1.1 求导数(差分)
  (1)3点公式
  (2)误差更大些的“2点公式”:前向差分公式,后向差分公式
  (3)5点公式
  下面是几个公式的比较,书中有77版代码。

点击(此处)折叠或打开

  1. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  2. ! chap1a.f90
  3. ! Computational Physics: Fortran Version
  4. ! 3-point formula,2-point formulas,5-point formula
  5. program main
  6. !ccccccccccccccccccccccccccccccccccccccccccccccccc!
  7. !
  8. !variables:
  9.  implicit none
  10.  integer :: counter ! using in do loop
  11.  integer,parameter :: num=100   ! using in do loop
  12.  real :: x=1.0,exact ! x: variable. exact: cos(1.0).
  13.  real :: threeP,twoFP,twoBP,fiveP   ! threeP: three-point. twoFP: two-point forward.
  14. ! twoBP: two-point backward. fiveP: five-point
  15.  real :: diff1,diff2,diff3,diff4 ! diff=cos(1.0)-approximate
  16.  real :: h=0.1 ! the step in calculation
  17. !
  18.  exact=cos(x)
  19. !
  20.  do counter=1,num,1
  21.   threeP=(sin(x+h)-sin(x-h))/(2*h)
  22.   diff1=exact-threeP
  23.   write(*,*)"exact:",exact,"threeP:",threeP,"diff1:",diff1
  24. !
  25.   twoFP=(sin(x+h)-sin(x))/h
  26.   diff2=exact-twoFP
  27.   write(*,*)"exact:",exact,"twoFP:",twoFP,"diff2:",diff2
  28. !
  29.   twoBP=(sin(x)-sin(x-h))/h
  30.   diff3=exact-twoBP
  31.   write(*,*)"exact:",exact,"twoBP:",twoBP,"diff3:",diff3
  32. !
  33.   fiveP=(sin(x-2*h)-8*sin(x-h)+8*sin(x+h)-sin(x+2*h))/(12*h)
  34.   diff4=exact-fiveP
  35.   write(*,*)"exact:",exact,"fiveP:",fiveP,"diff4:",diff4
  36.   h=h/2.0
  37.   write(*,*)"*****************","over of the",counter,"th line","**********************************"
  38. !
  39.  end do
  40. !
  41. stop
  42. 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行,误差才达到以上的程度。
  关于舍入误差以及其他误差来源,以后再专门学习一下吧。

阅读(889) | 评论(0) | 转发(0) |
给主人留下些什么吧!~~