FORTRAN数值方法及其在物理学中应用.ppt

数值方法及其在物理学中应用 School of Science, Xidian UniversitySchool of Science, Xidian University 1 课课 程程 提提 纲纲 第一章 算法语言简介与误差分析初步 1-1 FORTRAN 语言简介 1-2 MATLAB语言简介 1-3 简单力学问题的计算 1-4 误差及减小误差的原则 第二章 物理图形、图像与计算机模拟 2-1 简谐振动曲线及其合成 2-2 阻尼运动与阻尼振动的模拟 2-3 驻波的模拟 2-4 点电荷与点电荷系的等势线模拟 2-5 波的干涉和衍射图形模拟 2 第三章 物理学中定积分的数值计算方法 3-1 定积分基本数值算法及其应用 3-2 龙贝格法及其应用 3-3 高斯求积法 第四章 物理学中常微分方程初值问题的数值解法 4-1 物理学中的常微分方程 4-2 常微分方程初值问题的一级、二级欧拉近似法 4-3 龙格库塔法 第五章 物理学中线性代数方程组的数值解法 5-1 物理问题与线性代数方程组 5-2 高斯列主元消去法 5-3 解三对角方程组的追赶法 5-4 线性方程组的迭代解法 5-5 积分方程的数值解法 3 第六章 物理学中非线性方程的求根问题 6-1 根的搜索和二分法 6-2 函数迭代法 6-3 牛顿迭代法 第七章 实验物理学中的插值和数据拟合 7-1 实验数据的拉格朗日插值法 7-2 差商与牛顿插值公式 7.3 Hermite插值 7-4 三次样条插值 7-5 数值微分 7-6 最小二乘曲线拟合法 4 课课 程程 特特 点点 数值方法计算机编程大学物理知识数值方法计算机编程大学物理知识 注重利用有关数学结论求解物理问题结果，对注重利用有关数学结论求解物理问题结果，对 有关定理证明不作深入研究有关定理证明不作深入研究 5 课时：课时：4646学时学时2424机时机时 参考书：参考书： 1 施吉林等编. 计算机数值方法. 北京：高等教育出版社，2001。 2 王世儒等编. 计算方法. 西安：西安电子科技大学出版社，1996。 3 电子科技大学应用数学系编. 实用数值计算方法. 北京：高等教育出版社，2001。 4J.Stoer, R.Bulirsch. Introduction to Numerical Analysis. New York: Springer- Verlag, 1980. 5 李有法编. 数值计算方法. 北京：高等教育出版社，1996。 6 谭浩强，田淑清编著. FORTRAN语言. 北京：清华大学出版社，1990。 7 徐士良编著. FORTRAN常用算法程序集. 北京：科学出版社，1995。 8 刘云鹏编. 计算物理. 西安：西安电子科技大学教材科，1986。 9 吴百诗编. 大学物理. 西安：西安交通大学出版社，1994。 6 第一章 算法语言简介与误差分析初步 1.1 FORTRAN语言简介 1.2 MATLAB 语言简介 1.3 简单力学问题的计算 7 一、FORTRAN语言的三点说明 1.1 FORTRAN语言简介 1. 常数与变量的说明 2. 语句书写格式 Formula Translator(公式翻译器，数学语言翻译为机器语言） 3. 与C语言的不同 8 1. 常数与变量的说明 FORTRAN中常数与变量分为7类 n 实型（real）32bit； n 整数型（integer）长整型32bit，短整型16bit； n 双精度实型（double precision）real*8； 8bytes 64bits n 复型（complex）a+bi 两个浮点数；也分单双精度 n 双精度复型（D-P-C）complex*16； n 逻辑型（logical） True, False 分别以1，0代表； n 文字型（character） （字符型） 9 n 在FORTRAN 77中有字符型变量，字符型常数 关于变量的说明： n 定义注意I-N法则 n 名称长度不超过6个字符长（字母开头） n 大小写等价 n implicit语句： implicit real*8 (a, c), (t-v) implicit integer (d, e) 只能存储在字符型变量中。 10 2. 语句书写的格式 FORTRAN77扩展名为: *.f或*.for FORTRAN90无严格限制（以.f90为后缀） n 第1列有字符C时（注释行），不参加编译和运行。 FORTRAN90普遍用！放在语句后，但该句参与 本行在它后面的所有均被定义为注释。 n 当第6列上有非零和非空格的字符时，表明此行 为上行的继续行。例如：& n语句在772列中；语句标号15列。 运行。 11 3. 与C语言的不同 n 不分大小写。 n 每句末尾不必要写分号。 n 程序代码命令间的空格没有意义。 n 与C不同，FORTRAN不使用 。 n 数据类型多出了复数、乘幂运算和逻辑判断类型。 a=cmplx(1.0, 2.0) a=(1.0, 2.0) ! 1+2i b=3.0*(1./5.) ! 12 二、FORTRAN基本语句 1、可执行语句 (1) 赋值语句 V(变量)=e(表达式) (2) 流程控制语句 (3) 输入、输出语句 read(u,*), write(u,*) 13 ncall s(d1,d2,dn) （2）流程控制语句 n 无条件 goto语句 goto k n 算术条件语句 if(e) k1, k2, k3 (数轴记忆法) (e0) n 逻辑条件语句 if(e) s n 循环do语句 n 继续语句 continue return (在end前) do n i=m1,m2,m3 do n i=m1,m2 逻辑表达式 算术表达式 14 算数条件语句: if(e) k1, k2, k3 (e0) 例1：编程给出的值。 为键盘输入。要求： 15 算数条件语句的计算编程 cha1-1.f: read(*,40) X 40 format(F8.2) if(x)10,20,30 10 Y=-1.57079 goto 100 20 Y=0 goto 100 30 Y=1.57079 goto 100 100 write(*,50)x,Y 50 format(1X,2HX=,F10.6,4H, Y=F10.6) end 16 关系运算符号： FORTRAN 77: .gt. .ge. .lt. .le. .eq. .ne. FORTRAN 90 : = = = /= 逻辑条件语句 逻辑运算符号： .and. .or. .not.(逻辑非) .euv.(逻辑等) .neqv.(逻辑不等) if(X.le.2.1)Y=0.5*X+0.95 Y=0.7*X+0.53 write(*,*)X,Y 错误表示! 例2： 17 正确表示: if(X.le.2.1) Y=0.5*X+0.95 if(X.gt.2.1) Y=0.7*X+0.53 write(*,*)X,Y 或者 if(X.le.2.1) then Y=0.5*X+0.95 else Y=0.7*X+0.53 end if write(*,*)X,Y 18 cha1-3.f: read(*,20)x,y,z 20 format(3F10.4) big=x if(y.gt.big)big=y if(z.gt.big)big=z write(*,*)big=,big end 例3：有三个数x,y,z，要求打印出其中最大的数。 19 循环do语句 do n i=m1，m2，m3 或 do n i=m1，m2 即: do 标号 循环变量=表达式1，表达式2 ，表达式3 循环变量的初值 循环的终止值 循环增量值(省略则默认为1) 使用do语句循环时，可省略标号 do i=m1，m2，m3 循环程序 end do 例如： do 10 i=1,100,2 do 100 j=1, 8 do 60 x=1.2, 3.6, 0.2 do i=1，100 i=i+1 end do 例： 20 例4：编程求解0.0, 0.1, 0.2, 0.3的平方根 n do 10 i=0, 0.3, 0.1 x=0.0 do 10 i=1,4 y=sqrt(x) write(*,20)i,x,y 10 x=x+0.1 20 format(1X,I5,2F10.4) end 循环do语句 （错误） 或者用以下程序: do 10 x=0.0,0.3,0.1 y=sqrt(x) 10 write(*,20)i,x,y 20 format(1X, I5, 2F10.4) end 21 n call s(d1,d2,dn) return (在end前) n 继续语句 continue 22 (3) 输入、输出语句 read(u1,n1), write(u2,n2) 若u1，u2为“*”，则分别表示按键盘格式输入、屏幕表列 n表示格式说明: u1: 输入设备通道号:(17) u2: 输出设备通道号:(17) 输出。 若n=n1，则表示输入、输出按n1语句标号规定格式执行。 若n为“*”，则表示按自由格式输入或输出。 23 2. 非执行语句 (1)说明语句: 类型说明语句 、 隐含说明语句Implicit 例如 implicit integer(A,B,C) A,B,C开头的变量都是为整型 implicit none 关闭默认类型功能，所有变量都要事先声名 维数语句： dimension r(i)或r(i, j, k) real*8 a(i), complex*16 r(i) 公用语句： common w1, w2 主程序与子程序w1,w2同变量 参数语句： parameter v1=c, v2=c 24 (3)data语句(数据初值语句) n data v1/d1/,v2/d2/,vn/dn/ n或 data v1,v2,vn/d1,d2,dn/ (2)format语句 子例程子程序 subroutine s(a1,a2,an) (4)定义函数语句：f(a1,a2,an)=e (5)子程序语句: 函数子程序 function f(a1,a2,an) 25 函数子程序 function举例: integer fac,p,r write(*,*)n=,r=? read(*,*)n,r p=fac(n)/fac(n-r) write(*,*)n,r,p end integer function fac(n) fac=1 if(n.le.1) goto 77 do 10 k=2,n 10 fac=fac*k 77 return end 例5：编程求的值。 26 integer p,r write(*,*)n=,r=? read(*,*)n,r call fac(n,m) m1=m call fac(n-r,m) m2=m p=m1/m2 write(*,*)n,r,p end subroutine fac(i,m) m=1 if(i.le.1) goto 77 do 10 k=2,i 10 m=m*k 77 return end 子例程子程序 subroutine举例: 例5：编程求的值。 27 open(1,file=N!.dat) write(*,*)input N=? read(*,*)N M=1 I=2 5 M=M*I I=I+1 if(I.gt.N) goto 10 goto 5 10 write(1,*) M end 若不采用循环，还可以采用以下计算程序: 关于： 28 三、源程序语句排列顺序 n（1）说明语句（类型语句、维数语句等） n（2）数据语句（数据初值语句，定义函数语句） n（3）可执行语句 n（4）结束语句（格式语句可放在任意处） n（5）子程序语句（函数子程序，子例程子程序） 29 Aw以w个字符宽宽度输输出字符串 BN定义义文本框中的空位为为没有东东西，在输输入时时才需要使用 BZ定义义文本框中的空位代表0，在输输入时时才需要使用 Dw.d以w个字符宽宽来输输出指数类类型的浮点数，小数部分占d个字符宽宽 EW.dEe以w个字符宽宽来输输出指数类类型的浮点数，小数部分占d个字符宽宽 ，指数部分占e个字符 ENW.dEe 以指数类类型来输输出浮点数 ESW.dEe 以指数类类型来输输出浮点数 Fw.d以w个字符宽宽来输输出浮点数，小数部分占d个字符宽宽 Gw.dEe以w个字符宽宽来输输出任何种类类的数据 Iw.m以w个字符宽宽来输输出整数，最少输输出m个数字 Lw以w个字符宽宽来输输出T或F的真假值值 nX把输输出的位置向右跳过过n个位置 Format 格式 30 /代表换换行 ：在没有更多数据时结时结 束输输出 kPK值值控制输输入输输出的Scale Tn输输出的位置移动动到本行第n列 TLn输输出的位置向左相对对移动动n列 TRn输输出的位置向右相对对移动动n列 SP在数值为值为 正时时加上“正号” SS取消SP Bw.m把整数转换转换 成二进进制来输输出，输输出会占w个字符宽宽，固定输输 出m个数字， M值值可以不给给定 Ow.m把整数转换转换 成八进进制来输输出，输输出会占w个字符宽宽，固定输输 出m个数字， M值值可以不给给定 Zw.m把整数转换转换 成十六进进制来输输出，输输出会占w个字符宽宽，固定 输输出m个数字，M值值可以不给给定 31 四、FORTRAN常用内部函数（库函数） 功能 通用名 含义义功能通用名含义义 转换转换 到整 型 int(x) 正弦sin(x), dsin(x), csin(x) sin(x) 转换转换 到实实 型 float(x) 余弦cos(x), dcos(x), ccos(x) cos(x) 转换转换 到复 型 cmplx(a,b) a+ib 正切tan(x), dtan(x), ctan(x) tg(x) 平方根sqrt(x), dsqrt(x), cqrt(x) x1/2 反正弦asin(x), dasin(x)sin-1(x) 取绝对值绝对值abs(x), dabs(x), cabs(x) |x|反余弦acos(x), dcos(x) cos-1(x) 指数exp(x), dexp(x), cexp(x) ex 反正切atan(x), dtan(x) tg-1(x) 自然对对数alog(x), dlog(x) ,clog(x) lnx求共轭轭conjg(x) 常用对对数alog10(x), dlog(x) log10x 求余mod(a1,a2)a1/a2余数 最大值值max(a1,a2,a3) 最小值值min(a1,a2,a3) 32 错误 正确 n a(-b) a*-b a*(-b) n (ab)3 a*b*3 (a*b)*3 n sin6t sin6t sin(6*t) n aex a*e*x a*exp(x) n6.8log26.5 6.8*log26.5 6.8*alog10(26.5) 关系运算符号： +，*，/，* 表达式写法正误比对举例： n 33 1. “/”及“*”不能省略；一律用小括号。 2同类型算术量间才能运算。 说明 3括号 函数 * *，/ +， 34 五、有关循环语句 n（1）循环变量最好不要第二次赋值。 n（2）可以从循环体内转向循环体外， 不允许循环体外转向循环体内。 n（3）注意多重循环问题。 35 implicit real*8(A-H,O-Z) open(1,file=sin.dat) write(*,*)input w,N read(*,*)w,N pi=3.1415926 do 10 I=1,N t=2.*pi/w t=t*float(I)/N x=sin(w*t) write(*,*)t, x 10 write(1,*)t, x end 例6：编程给出的数据文件。 36 n（指数增加每次加1，可以用i代替） a=1 do 10 i=1, 63 10 a=a+2*i write(*,*) a end 例7：编程计算： a=1 do i=1, 63 a=a+2*i enddo write(*,*) a end 37 write(*,*)input x, N read(*,*)x,N y=0 do 20 M=1,N t=1 s=x do 100 J=1,M t=t*float(J) 100 s=s*x y=y+s/t 20 write(*,*) M,y end 例8：编程计算 的值。 y=0 do M=1,N t=1 s=x do J=1,M t=t*float(J) s=s*x enddo y=y+s/t write(*,*) M,y enddo end 38 六、FORTRAN语言的优缺点 n历史长，正确可靠，兼容性强。 n语法书写要求严格，更适合严谨的科学计算领域。 n可以直接对数组和复数进行运算。 n在并行方面，Fortran语言仍是不可替代的。 nFortran语言是一种编译语言，运行速度快。 nFortran语言自身仍在不断完善和发展，功能不断增强。 优点： 39 v本课程后续章节中，有关例题均采用Fortran编译， 如若用其他算法语言编程，Fortran源码仅供参考。 缺点： FORTRAN语言是编译语言，具有执行速度快等特点。 如格式严格必然导致灵活性下降，操作指针困难等等。 也正是这些优点带来了一些缺点。 40 EX1-1: 编程计算 ，直至 大于 为止。 EX1-2: 编程求 值， 由键盘输入 EX1-3: 编程计算 的和，直 至最后一项小于 。 作业 41 1.2 MATLAB 语言简介 n观看教学录像 片 42 一、瞬时性与极限 例9：已知某一质点沿某一直线运动，位移与时间满足方程 ： 求: t=2时,=1, 0.1, 0.01, 0.001, 0.0001的 及 时的瞬时值。 1.3 简单力学问题的计算 瞬时值： 43 read(*,*) t do 10 K=1,5 dt=10.*(1-K) v1=9.*t+4.5*dt-6.*t*t-6.*t*dt-2.*dt*dt a1=9.-12.*t-6.*dt 10 write(*,*)dt,v1,a1 v=9.*t-6.*t*t a=9.-12.*t write(*,*)v,a end 计算程序 44 二、运动方程问题 n 解: 例10：已知时有 。 求：时的 。 45 read(*,*) ts, N dt=ts/float(N) t=0.0 x1=0.0 do 10 I=1,N v=5.-3.*t*t x1=x1+v*dt 10 t=t+dt write(*,*)N,x1 end 输入：1.0, 500 运行结果：N=500, x1=4.003006 计算程序 46 弦截法 本问题等价于（求 这一超越方程的解） 弦截法求方程 根问题： （语句 ） 例11：已知求：时的 47 弦截法 直至（如 ） 用 代替两点中与其同号的点，得到一个缩小的区间 。 有取， 48 弦截法 （语句 ） 49 弦截法 （语句 ） 50 （语句 ） 弦截法 51 read(*,*)a,b,v,e fa=a+alog(1.+a)-v fb=b+alog(1.+b)-v if(fa*fb.gt.0) goto 1000 if(abs(fa).gt.e) goto 5 c=a goto 100 5 if(abs(fb).gt.e) goto 10 c=b goto 100 10 c=b-(a-b)/(fa-fb)*fb fc=c+alog(1.+c)-v if(abs(fc).le.e) goto 100 if(fa*fc.gt.0) goto 15 fb=fc b=c !c与b同号 goto 10 15 fa=fc a=c !c与a同号 goto 10 100 write(*,*) c=, c goto 111 1000 write(*,*) No solution 111 stop end 输入：0.0,1.0, 1.0, 1E-6 运行结果：t=5.571457E-1 计算程序 考虑程序用数组及循环？ 52 EX1-4: 质量为m的快艇，以 的速率行驶，发动机关闭后， 受到的阻力为 ，已知 m/s, m-1， 求15秒后快艇的速率和所走的路程各为多少？ EX1-5: 质点沿半径为0.1m的圆周运动，其角位置 表示为: (rad) 编程计算何时切向加速度和法向加速度恰好有相等的值 (误差小于 )。 作业 53 一、误差及其分类 1.4 误差及减小误差的原则 1. 模型误差：由实际问题抽象、简化为数学问题 2. 观测误差：测量工具限制或在数据的获取时随机因素 因此数值方法中的取数和运算往往会出现误差。 （建立数学模型时）所引起的误差。 所引起的物理量的误差。 在计算机上参与运算的数必须是有限小数或整数。 54 3. 截断误差：是用数值法求解数学模型时得到的正确解和 例如若用级数： 截断误差为： 4. 舍入误差：由于计算机所表示的位数有限，通常用 模型准确解间的误差方法误差。 的近似值，的前三项计算即取 四舍五入的办法取值而引起的误差。 55 数值计算方法所考虑的主要是计算误差。 56 二、绝对误差和相对误差 绝对误差： （ 为近似值， 为准确值） 相对误差： 例：某甲测量一物体长度为 cm， cm， 则有： 某乙测量一物体长度为 cm， cm 57 解：球体积积的计计算公式为为： 当半径有一误差时，球体积的相对误差为： 现要求， 例：计计算球体积积要使相对误对误 差限为为1%,问问度量半径 时允许的 相对误差限是多少？ 即 58 三、有效数字 定义：若近似值 的绝对误差限是某一位上的半个单位 ，该位到 的非零数字一共有 位，则称近似值 具 有 位有效数字。 按四舍五入取三位小数得的近似值为3.14，取六位小数 v 对于同一个数的近似值而言，有效数字位数越多， 例：3.14159265 的近似值为3.14159。无论取几位小数所得近似值，其绝对 误差都不会超过其末位数的半个单位。 其绝对误差与相对误差都越小； v 绝对误差或相对误差越小，有效数字的位数有可能越多。 59 四、数值计算中减小误差的原则 1选用数值稳定性好的算法 例12：考虑循环公式 随着 的增加，误差迅速增加。 设 的误差为 ，求 值时的误差为多少？ 考察引起的误差是否累计增长。如不增长就认为是稳 定的；严重增长就认为不稳定。 60 “蝴蝶效应” 蝴蝶效应:气象学家洛伦兹1963年提出。其大意为： 其原因在于：蝴蝶翅膀的运动，导致其身边的空气 系统发生变化，并引起微弱气流的产生，而微弱气 流的产生又会引起它四周空气或其他系统产生相应 的变化，由此引起连锁反映，最终导致其他系统的 极大变化。 偶尔扇动几下翅膀，可能在两周后引起美国德克萨 斯引起一场龙卷风。 一只南美洲亚马孙河流域热带雨林