物理类模拟

XJTU 物理类模拟 2009022023 王梦冰 2011/8/3 1 / 10 目录 工贸企业有限空间作业目录特种设备作业人员作业种类与目录特种设备作业人员目录1类医疗器械目录高值医用耗材参考目录 ..................... Contents 用 Fortran 模拟单摆运动 .............................. 2 实验原理 ...................................................... 2 实验过程 ...................................................... 4 附：程序源代码： ...................................... 4 分形问题 ........................................................ 7 Mathematica程序源代码： ....................... 7 光子散射后的能量抽样问题： .................. 10 2 / 10 用 Fortran 模拟单摆运动 实验原理 解二阶微分方程的 Runge-Kutta法： )',,('' yyxfy  ，此处 )',,( yyxf 是子程序。 Function ydpri (xn,yn,dyn,iflag) Xn 为自变量，yn 为因变量，dyn 是一阶倒数，iflag 为通信标志。 Desork(iter,h,x,y,dy,iflag) 形参名 用途 说明 Iter 输入 带球的个数 h 输入 步长，即 x(n)和 x(n+1)的间隔 x 输入 自变量 y 输出 因变量 dy 输出 y 的一阶导数 iflag 输入 在函数子程序 ydopri 中指明 要求解哪一位的通信标志 设 ),,( uyxf dx dy  ),,( uyxg dx du  ),,(0 unn uyxhfk  ),,(0 nnn uyxhgm  ) 2 , 2 , 2 ( 001 m u k y h xhfk nnn  ) 2 , 2 , 2 ( 001 m u k y h xhgm nnn  ) 2 , 2 , 2 ( 112 m u k y h xhfk nnn  3 / 10 ) 2 , 2 , 2 ( 112 m u k y h xhgm nnn  ),,( 223 mukyhxhfk nnn  则：  32101 22 6 1 kkkkyy nn   32101 22 6 1 mmmmuu nn  设  yxf dx dy yu ,'   ',,'' yyxgy dx du  带入上面的公式，得： nhyk '0   nnn uyxhgm 0        2 0' 1 m yhk n ) 2 , 2 , 2 ( 001 m u k y h xhgm nnn         2 1' 2 m yhk n ) 2 , 2 , 2 ( 112 m u k y h xhgm nnn   2''3 myhk  ),,( 223 mukyhxhgm nnn   32101 22 6 1 kkkkyy nn  4 / 10  3210'' 1 22 6 1 mmmmyy nn  实验过程 经子程序和主程序，得： 当角度接近 0 时， 430175724.00  E 4900632586.00  E T/4=   2 )()( 00    tt =12.1061082s T=4*T/4=48.4244328 附：程序源代码： subroutine desork(iter,h,x,y,dy,iflag)!解二阶微分方程的子程序 C subroutine to solve second-order ordinary differential C equations y''=f(x,y,y') by runge-kutta method. dimension x(121080),y(121080),dy(121080) do 10 i=1,iter-1 tko=h*dy(i) tmo=h*ydopri(x(i),y(i),dy(i),iflag) xn=x(i)+h/2. yn=y(i)+tko/2. dyn=dy(i)+tmo/2. tk1=h*(dy(i)+tmo/2.) tm1=h*ydopri(xn,yn,dyn,iflag) yn=y(i)+tk1/2. dyn=dy(i)+tm1/2. tk2=h*(dy(i)+tm1/2.) tm2=h*ydopri(xn,yn,dyn,iflag) xn=x(i)+h yn=y(i)+tk2 dyn=dy(i)+tm2 tk3=h*(dy(i)+tm2) tm3=h*ydopri(xn,yn,dyn,iflag) y(i+1)=y(i)+(tko+2.*tk1+2.*tk2+tk3)/6. dy(i+1)=dy(i)+(tmo+2.*tm1+2.*tm2+tm3)/6. 10 continue end 5 / 10 C program solution ! 主程序 dimension x(121080),y(121080),dy(121080)!定义大小为 121080 的数组 real x,y,dy!限定变量都为实数 integer::n n=2 x(1)=0. do while(n<=121080) x(n)=x(n-1)+0.0001 n=n+1 end do!用 do 循环输入自变量 x iter=121080 h=0.0001 iflag=1!从 x 数组第一个元素开始算起 y(1)=30.!y 的初值 dy(1)=0.!dy 的初值 call desork(iter,h,x,y,dy,iflag) do 10 i=1,iter write(*,5) x(i),y(i) 10 continue 5 format(1x,2e18.9)!依次输出 x,y end C function ydopri(xn,yn,dyn,iflag)!所要解的微分方程 ydopri=-sind(yn) end 6 / 10 7 / 10 分形问题 Mathematica 程序源代码： redokoch[ptlist_List]:= Block[{tmp={},i,pnum=Length[ptlist]}, For[i=1,i 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 为：                                22 2 11111 )( 1 )( xxxx x K xf     1 2 2 2 1 2 2 1 2 1 ( ) 2 1 ( ) 2 2 1 1 ( ) ( )(1 2 ) 2 1 ( ) ( ) f x x f x x H x K x x H x K x                                  )()()()()( 2211  xfxHxfxHxf  1 2 2 2 1 1 2 2 2 2 1 1 1 2 1 2 1 2 1 2 4 ( ) ( )(1 2 ) 8 ( ) 27 ( ) 2 (1 2 ) 2(1 )(1 2 ) ln(1 2 ) 2 (1 ) ( ) ( ) (1 2 ) f f X X H x M K H x M K P H x f x dx                                        2 2 2 2 2 )( ff f XX M XH      > 1 1 1 1 2 )( ff f XX M XH      11 P >