南开大学计算物理讲义7

1、In this section we look at how to solve higher order ordinary differential equations. Some of the techniques discussed here are based on the techniques introduced in 7, and some rely on a combination of this material and the linear algebra of 5.8 Higher order ordinary differential equations8.2.1 Ini

2、tial value problems (初始值问题初始值问题)8 Higher order ordinary differential equations8.2.2 Boundary value problems (边值问题边值问题) 8.1 First order ordinary differential equation sets 8.2 Higher order ordinary differential equations(1) The Taylor Series expansion0)0(, 1)0(423232yxttyxytttyxx ! 3 !2)( ! 3 !2)(323

3、2yhyhhyyhtyxhxhhxxhtx8.1 First order ordinary differential equation sets 2322 ttyxx323242ttyxytttyxxtyxx62 238 ttyxytyxy68 6 6 yxyyxx ! 3 !2)( ! 3 !2)(3232yhyhhyyhtyxhxhhxxhtx6 6 yxyyxxh=0.01, T=0 2238 322 ttyxyttyxx323242ttyxytttyxx6 6 yxyyxx ! 3 !2)( ! 3 !2)(3232yhyhhyyhtyxhxhhxxhtxX1=X-Y+T*(2.0-T

4、*(1.0+T)Y1=X+Y+T*T*(4.0+T)X2=X1-Y1+2.0-T*(2.0+3.0*T)Y2=X1+Y1+T*(8.0+3.0*T)X3=X2-Y2-2.0-6.0*TY3=X2+Y2-8.0+6.0*TX4=X3-Y3-6.0Y4=X3+Y3+6.0X=X+H*X1+H*(X2/2.+H*(X3/6.+H*X4/24.)Y=Y+H*Y1+H*(Y2/2.+H*(Y3/6.+H*Y4/24.)T=REAL(K)*Htyxytyxx68 62 0)0(, 1)0(yxDATA T/0.0/,X/1.0/,Y/0.0/,H/0.01/,NSTEP/100/DO K=1,NSTEP2

5、238 322 ttyxyttyxx323242ttyxytttyxx6 6 yxyyxx ! 3 !2)( ! 3 !2)(3232yhyhhyyhtyxhxhhxxhtxX1=X-Y+T*(2.0-T*(1.0+T)Y1=X+Y+T*T*(4.0+T)X2=X1-Y1+2.0-T*(2.0+3.0*T)Y2=X1+Y1+T*(8.0+3.0*T)X3=X2-Y2-2.0-6.0*TY3=X2+Y2-8.0+6.0*TX4=X3-Y3-6.0Y4=X3+Y3+6.0X=X+H*X1+H*(X2/2.+H*(X3/6.+H*X4/24.)Y=Y+H*Y1+H*(Y2/2.+H*(Y3/6.+H

6、*Y4/24.)T=REAL(K)*HWRITE(*,*)T,X,YENDDOSTOPENDtyxytyxx68 62 0)0(, 1)0(yx0)0(, 1)0(423232yxttyxytttyxx32)sin()()cos()(ttetyttetxtt解析解解析解: Euler methodTnnTnnxxxxXxxxxX, , , ,121121)(XFX ),(),(),(2121222111txxxfdtdxtxxxfdtdxtxxxfdtdxnnnnn)Y,t (YY1nnnnhf)X,t (X1nnnnhfX x1221kxxxxTnnTnnxxxxXxxxxX, , , ,1

7、21121)(XFX )X,t (X1nnnnhfX 1221kxxxx)(XFX )X,t (X1nnnnhfX 11 . 0011 . 010X199. 01 . 01 . 011 . 011 . 0X2例例98. 0199. 01 . 099. 01 . 099. 01 . 0X3323242ttyxytttyxx 3121323312113221421xxxxxxxxxxxx x1=t, x2=x, x3=y3121323121132421xxxxxxxxxX)(XFX 例例Yn+1 = Yn + (k(1) + 2k(2) + 2k(3) + k(4)/6.Yn+1 = Yn + (

8、k(1) + k(2)/2.k(2)= h f(tn +h, Yn + k(1)k(1)=h f(tn, Yn);k(1)=h f(tn, Yn);k(2)= h f(tn +h/2, Yn + k(1)/2)k(3)= h f(tn +h/2, Yn + k(2)/2)k(4)= h f(tn +h, Yn + k(3)Runge-Kutta methodsXn+1 = Xn + h(F(1) + 2F(2) + 2F(3) + F(4)/6.F(1)= F(Xn);F(2)= F( Xn + hF(1)/2)F(3)= F(Xn + hF(2)/2)F(4)= F(Xn + hF(3)(X

9、FX Yn+1 = Yn + h(k(1) + 2k(2) + 2k(3) + k(4)/6.k(1)= f(tn, Yn);k(2)= f(tn +h/2, Yn +h k(1)/2)k(3)= f(tn +h/2, Yn + hk(2)/2) k(4)= f(tn +h, Yn + hk(3)0)0(, 1)0(423232yxttyxytttyxx; 0) 0(; 1) 0(; 0) 0(4213213121323312113221xxxxxxxxxxxxxxxx1=t, x2=x, x3=y例例F(1)= F(Xn);F(2)= F( Xn + hF(1)/2)F(3)= F(Xn +

10、 hF(2)/2)F(4)= F(Xn + hF(3)Xn+1 = Xn + h(F(1) + 2F(2) + 2F(3) + F(4)/6.; 0) 0(; 1) 0(; 0) 0(4213213121323312113221xxxxxxxxxxxxxxxx1=t, x2=x, x3=yF(1)= F(Xn);F(2)= F( Xn + hF(1)/2)2111401112011) 1 (F h=0.01 Xn + hF(1)/2=01. 0005. 1005. 02/ ) 2(*01. 002/01. 012/01. 0; 0) 0(; 1) 0(; 0) 0(42132131213233

11、12113221xxxxxxxxxxxxxxxF(1)= F(Xn);F(2)= F( Xn + hF(1)/2)3055. 02005. 0*401. 0005. 13005. 02005. 0005. 0*201. 0005. 11)2(F Xn + hF(1)/2=01. 0005. 1005. 02/ ) 2(*01. 002/01. 012/01. 0F(1)= F(Xn);F(2)= F( Xn + hF(1)/2)F(3)= F(Xn + hF(2)/2)F(4)= F(Xn + hF(3)Xn+1 = Xn + h(F(1) + 2F(2) + 2F(3) + F(4)/6.可

12、使用与求可使用与求 F(2)相同的方法求相同的方法求 F(3) 和和 F(4)SUBROUTINE XPSYS ( X, F); 0) 0(; 1) 0(; 0) 0(4213213121323312113221xxxxxxxxxxxxxxxDIMENSION X(3), F(3) F(1)=1.0F(2)=X(2)-X(3)+X(1)*(2.0-X(1)*(1.0+ X(1)F(3)=X(2)+X(3)-X(1)* X(1)*(4.0-X(1)RETURNENDSUBROUTINE RK4SYS(N, T, X, H, NSTEP)DIMENSION X(10),Y(10), F1(10),

13、 F2(10), F3(10), F4(10)DO K=1, NSTEP CALL XPSYS(X, F1) DO I=1, N Y(I)=X(I)+H*F1(I)/2. ENDDOENDDO RETURN CALL XPSYS(Y, F2) DO I=1, N Y(I)=X(I)+H*F2(I)/2. ENDDO CALL XPSYS(Y, F3) DO I=1, N Y(I)=X(I)+H*F3(I) ENDDO CALL XPSYS(Y, F4)DO I=1, N X(I)=X(I)+H*(F1(1)+2.*(F2(I)+F3(I)+F4(I)/6. ENDDOT=T+REAL(K)*H

14、WRITE(*,*) T, (X(I), I=1,N)dimension x(3)data n/3/,t/0.0/,x/0.0,1.0,0.0/,h/0.01/,nstep/100/call rk4sys(n,t,x,h,nstep)stopend; 0) 0(; 1) 0(; 0) 0(4213213121323312113221xxxxxxxxxxxxxxx8.2.1 Initial value problems (初始值问题初始值问题)8.2 Higher order ordinary differential equationsy = f(t, y, y) subject to y(0

15、) = c0 and y(0) = c1 subject to y(0) = c0 and y(1) = c1 Initial value problems (初始值问题初始值问题) Boundary value problems (边值问题边值问题) y(1)=0, y(0)=0y(0)=0, y(1)=0where at t = t0 we know the values of x0, x1, ,xn-1 mtxFdtxd),(2210 xx mtxFx),(0110 xxxx13)0( , 7)0( , 3)0(sincos 2 yyyteyyyy13)0(, 7)0(, 3)0(sin

16、cos210210221102xxxtexxxxxxxx13)0(, 7)0(, 3)0(sincos210210212xxxtexxxxXx0021)(),(iiniiytyyyytfyx1=y1, x2=y2 , xn=yn)(XFX TnnTnnxxxxXxxxxX, , , ,121121TnnfxxxxXF,)(1210021)(),(iiniiytyyyytfyx1=t, x2=y1, x3=y2 , xn+1=yn)(XFX TnnTnnxxxxXxxxxX, , , , ,121121TnnfxxxxXF,)(121)(XFX where at t = t0 we know t

17、he values of y1, y2, , yn x1=t, x2=y1, x3=y2 , xn+1=yn0021)(),(iiniiytyyyytfyTnnxxxxX)0(),0(,),0(),0()0(12113)0( , 7)0( , 3)0() sin()cos( 2 yyyteyyyyx1=t, x2=y, x3=y , x4=y232434sincos1xexxxxXxTX13,7 , 3 ,0)0()/log()sin() cos() ( 2tytyyyyyHigher order equation sets 8.2.2 Boundary value problems (边值问

18、题边值问题) For second (and higher) order odes, two (or more) initial/boundary conditions are required. If these two conditions do not correspond to the same point in time/space, then the simple extension of the first order methods outlined in section 8.1 can not be applied without modification. 8.2 High

19、er order ordinary differential equationsThere are two relatively simple approaches to solve Boundary value problems.y = f(t, y, y) subject to y(0) = c0 and y(0) = c1 subject to y(0) = c0 and y(1) = c1 Initial value problems (初始值问题初始值问题) Boundary value problems (边值问题边值问题) y(1)=0, y(0)=0y(0)=0, y(1)=0

20、8.2.2.1 Shooting methody = f(t, y, y) With the shooting method we apply the y(0)=c0 boundary condition and make some guess that y(0) = a a0. This gives us two initial conditions so that we may apply the simple time-stepping methods already discussed in section 8.1. The calculation proceeds until we

21、have a value for y(1). subject to y(0) = c0 and y(1) = c1 If this does not satisfy y(1) = c1 to some acceptable tolerance, we revise our guess for y(0) to some value a a1, say, and repeat the time integration to obtain an new value for y(1). This process continues until we hit y(1)=c1 to the accepta

22、ble tolerance. The number of iterations which will need to be made in order to achieve an acceptable tolerance will depend on how good the refinement algorithm for a a is. 1 Newton-Raphson Setting the quadratic and higher terms to zero and solving the linear approximation of f(x) = 0 for x gives)(!)

23、()( ! 2)()( )()()()(02000 xfnxxxfxxxfxxxfxfnn0)( )()()(00 xfxxxfxf0)( )()()(11nnnnnxfxxxfxf)( / )(1nnnnxfxfxxWe may use the root finding methods discussed in section 3 to undertake this refinement. (chord) )( / )(1nnnnxfxfxx)/()(111nnnnnnnfffxxxx)(!)( !2)( )()()(2xfnxxfxxxfxfxxfnnxxfxxfxfkkxk)()(lim

24、)( 0)(kkxff )(1hOhfffkkkxhy= - 2 2 (y+1)/4(y+1)/4y(0)=0, y(0)=1y(0)=0, y(1)=1y(0)= y(1)=0, y(0)=0Here is a parameter to be adjusted in order to have|F( )|= |y (1) - 1| y(0)=0, y(1)=0 y (1) y(1)=1x3(0)= + Runge-Kutta methodsx1 =1x2 = x3x1(0)=0, x2(0)=0, x2(1)=1 x1=t, x2=y, x3=yy=- 2 (y+1)/4y(0)=0, y(

25、1)=1|F( )|= |x2 (1 )-1| x3 = -2 2( (x2+1)/4+1)/4x3(0)=( x2(1)- x2(0)/(1.-0.) )/()(111nnnnnnnfffxxxxX30=(YU-YL)/ (XU-XL)X31=X30+DX将将x1(0), x2(0) 和和x3(0)代入代入Runge-Kutta x2 (1 ) =f0将将x1(0), x2(0) 和和x3(0)+h代代入入Runge-Kutta x2 (1 ) =fnRETURNEND)/()(111nnnnnnnfffxxxx|F( )|= |x2 (1 )-1| dimension y0(3), y(3

26、)data n/3/,t/0.0/,y0/0.0,0.0,1.0/,nstep/100/ ,NMAX/10/ tol=1.0E-06 XL=0.0XU=1.0DX=0.01H=(XU-XL)/(NSTEP-1)YL=0.0YU=1.0|F( )|= |y 2 (1 ) - 1|X0=(YU-YL)/ (XU-XL)X1=X0+DXDO K=1, NMAXy=y0y(3)=X0F0=y(2)-1.0 y=y0 y(3)=X1F1=y(2)-1.0 D=F1-F0 STOPELSEENDIFENDDOENDcall rk4(n,t,y,h,nstep)call rk4(n,t,y,h,nstep)

27、SUBROUTINE XPSYS ( X, F)DIMENSION X(3), F(3) PI=4.0*ATAN(1.0)F(1)=1.0F(2)= X(3)F(3)=-PI*PI*(X(2)+1.0)/4.0RETURNENDx1 =1x2 = x3x3 = -2 2( (x2+1)/4+1)/4Y(X)=COS(x/2)+2x/2)+2SIN(x/2)-1x/2)-1y=- 2 2 (y+1)/4(y+1)/4y(0)=0, y(1)=1y = f(t, y, y) subject to y(0) = c0 and y(1) = c1 y(1)=0, y(0)=0y(0)=0, y(1)=

28、08.2.2.2 Linear equations The alternative is to rewrite the equations using a finite difference approximation with step size t = (t1t0)/N to produce a system of N+1 simultaneous equations. Consider the second order linear systemy + ay + by = c, with boundary conditions y(t0) = a a. and y(t1) + y(t1)

29、 = b b. yi (Yi+1 Yi-1)/2 t, yi (Yi+1 2Yi + Yi-1)/ t2 Substitution intoy + ay + by = c, If we use the central difference approximations t = (t1t0)/Ni = 1-N-1y + ay + by = c, (Yi+1 2Yi + Yi-1)/ t2+a(Yi+1 Yi-1)/2 t + b Yi= c Gives(1-1/2a t)Yi-1 + (b t22)Yi + (1+1/2a t)Yi+1 = c t2,(1-1/2a t)Y0 + (b t22)

30、Y1 + (1+1/2a t)Y2 = c t2,and we may write the boundary condition at Y0 = a a. y(t0) = a a and y(t1) + y(t1) = b b. Yn + Y n= b b , as t = t0,(1-1/2a t)Yi-1 + (b t22)Yi + (1+1/2a t)Yi+1 = c t2,(1-1/2a t)Y1 + (b t22)Y2 + (1+1/2a t)Y3 = c t2,(1-1/2a t)Yn2 + (b t22)Yn1 + (1+1/2a t)Yn = c t2,i = 1-N-1)(!

31、)( !2)( )()()(2xfnxxfxxxfxfxxfnn)(1hOhfffkkktFor t = t1 we must express the derivative in the boundary condition in finite difference form using points that exist in our mesh. The obvious thing to do is to take the backward difference so that the boundary condition becomesYn + (Yn Yn-1)/ t= b b , bu

32、t as we shall see, this choice is not ideal.y(t1) + y(t1) = b b.The above strategy leads to the system of equations Yn + Y n= b b , Y0 = a a,(1-1/2a t)Y0 + (b t 22)Y1 + (1+1/2a t)Y2 = c t 2,(1-1/2a t)Y1 + (b t 22)Y2 + (1+1/2a t)Y3 = c t 2,(1-1/2a t)Yn2 + (b t 22)Yn1 + (1+1/2a t)Yn = c t 2, Yn-1 + (

33、t + 1) Yn = bbtThis tridiagonal system may be readily solved using the method discussed in section 4.5. Tridiagonal matricesThere is one problem with the scheme outlined above: while the Y and Y derivatives within the domain are second order accurate, the boundary condition at t = t1 is imposed only

34、 as a first-order approximation to the derivative. yi (Yi+1 Yi-1)/2 t, yi (Yi+1 2Yi + Yi-1)/ t2 Yn + (Yn Yn-1)/ t= b b , y(t1) + y(t1) = b b.The first way is the most obvious: make use of the Yn-2 value to construct a second-order approximation to YtWe have two ways in which we may improve the accur

35、acy of this boundary condition.so that we have)( ! 3)( !2)( )()(32xfxxfxxxfxfxxf)( !3)2()( !2)2()( 2)()2(32xfxxfxxxfxfxxfYn ( 3Yn -4Yn-1 +Yn-2)/2 t, y(t1) + y(t1) = b b.Yn-2 4Yn-1 + (2 t + 3) Yn = 2bbtY0 = a a,(1-1/2a t)Y0 + (b t22)Y1 + (1+1/2a t)Y2 = c t2,(1-1/2a t)Y1 + (b t22)Y2 + (1+1/2a t)Y3 = c

36、 t2,(1-1/2a t)Yn2 + (b t22)Yn1 + (1+1/2a t)Yn = c t2,Yn-2 4Yn-1 + (2 t + 3) Yn = 2bbtThe problem with this is that we then loose our simple tridiagonal system.The second approach is to artificially increase the size of the domain by one mesh point and impose the boundary condition inside the domain.

37、 The additional mesh point is often referred to as a dummy(虚拟的虚拟的) mesh point or dummy cell. t)(!)( !2)( )()()(2xfnxxfxxxfxfxxfnn)(!)( !2)( )()()(2xfnxxfxxxfxfxxfnn)(2211hOhfffkkky(t1) + y(t1) = b b. Yn-1 + 2 t Yn +Yn+1= 2bbt tY0 = a a,(1-1/2a t)Y0 + (b t22)Y1 + (1+1/2a t)Y2 = c t2,(1-1/2a t)Y1 + (b

38、 t22)Y2 + (1+1/2a t)Y3 = c t2,(1-1/2a t)Yn2 + (b t22)Yn1 + (1+1/2a t)Yn = c t2,The attraction of this approach is that we maintain our simple tri-diagonal system.(1-1/2a t)Yn1 + (b t22)Yn + (1+1/2a t)Yn+1 = c t2, Yn-1 + 2 t Yn +Yn+1= 2bbt Clearly the same tricks can be applied if we have a boundary condition on the derivative at t0.Higher or