PROGRAM TESTAB C C---------------------------------------------------------------------- C C This is a test program for the two subroutines KUTSIM and C ADBASH, which may be used to find a solution to a system C of first order differential equations. KUTSIM uses the C Kutta-Simpson method and ADBASH uses the Adams-Bashford C method. When ADBASH is used, the first three steps are C calculated using KUTSIM. C C The system of equations for this test program is C dy(1)/dx = (a+b+1/x)y(1)/2 + (a-b+1/x)y(2)/2 C dy(2)/dx = (a-b+1/x)y(1)/2 + (a+b+1/x)y(2)/2 C which has the solution C y(1) = x exp(ax) + exp(bx) C y(2) = x exp(ax) - exp(bx) C C Written by Robert Walraven 9 Oct 81 C Applied Science C University of California, Davis 95616 C (916) 752-0360 C C---------------------------------------------------------------------- C DIMENSION YIP1(2),YI(2),YIM1(2),YIM2(2),YIM3(2),F(2),G(2) COMMON /CONST/ A,B LOGICAL ADBFLG C C Initialize variables, get information, and print titles C A = 2. B = 0.5 X = 0. YIM3(1) = 1. YIM3(2) = -1. NVARS = 2 Y1TRUE = YIM3(1) Y2TRUE = YIM3(2) WRITE (5,10) 10 FORMAT(' Test Program for subroutine KUTSIM and ADBASH'/ 1 ' Enter step size (You might try 0.01): '$) READ (5,*) H WRITE (5,20) 20 FORMAT(' Do you want to try Kutta-Simpson (K) or ', 1 'Adams-Bashford (A) ? '$) READ (5,30) IANS 30 FORMAT (A1) ADBFLG = IANS.EQ.'A' WRITE (5,40) 40 FORMAT(1X//9X'X'9X'Y1'12X'Y1,TRUE'11X'Y2'11X'Y2,TRUE'/) WRITE (5,50) X,YIM3(1),Y1TRUE,YIM3(2),Y2TRUE 50 FORMAT(1X,F10.3,4(1PE16.7)) C C Always do the first three steps by Kutta-Simpson method. C CALL KUTSIM (H,X,NVARS,YIM3,YIM2,F,G) X = X+H CALL YTRUE (X,Y1TRUE,Y2TRUE) WRITE (5,50) X,YIM2(1),Y1TRUE,YIM2(2),Y2TRUE CALL KUTSIM (H,X,NVARS,YIM2,YIM1,F,G) X = X+H CALL YTRUE (X,Y1TRUE,Y2TRUE) WRITE (5,50) X,YIM1(1),Y1TRUE,YIM1(2),Y2TRUE CALL KUTSIM (H,X,NVARS,YIM1,YI,F,G) X = X+H CALL YTRUE (X,Y1TRUE,Y2TRUE) WRITE (5,50) X,YI(1),Y1TRUE,YI(2),Y2TRUE C C Do the rest of the steps by the specified method. C DO 80 I=4,100 IF (ADBFLG) GO TO 60 CALL KUTSIM (H,X,NVARS,YI,YIP1,F,G) YI(1) = YIP1(1) YI(2) = YIP1(2) GO TO 70 60 CALL ADBASH (H,X,NVARS,YIP1,YI,YIM1,YIM2,YIM3,F,G) 70 X = X+H CALL YTRUE (X,Y1TRUE,Y2TRUE) WRITE (5,50) X,YIP1(1),Y1TRUE,YIP1(2),Y2TRUE 80 CONTINUE C C Th-th-th-that's all folks! C CALL EXIT END SUBROUTINE YTRUE (X,Y1,Y2) C C---------------------------------------------------------------------- C C This subroutine calculates the true values of y1 and y2 C for this test program. C C---------------------------------------------------------------------- C COMMON /CONST / A,B Y1 = X*EXP(A*X) + EXP(B*X) Y2 = X*EXP(A*X) - EXP(B*X) RETURN END SUBROUTINE FUNC (F, X, Y) C C---------------------------------------------------------------------- C C Test function subroutine for solution of differential C equations. The system to be solved is C dy(1)/dx = (a+b+1/x)y(1)/2 + (a-b+1/x)y(2)/2 C dy(2)/dx = (a-b+1/x)y(1)/2 + (a+b+1/x)y(2)/2 C which has the solution C y(1) = xexp(ax) + exp(bx) C y(2) = xexp(ax) - exp(bx) C C---------------------------------------------------------------------- C DIMENSION F(2), Y(2) COMMON /CONST/ A,B IF (X.EQ.0.) GO TO 10 T1 = 0.5*(A+B+1./X) T2 = 0.5*(A-B+1./X) F(1) = T1*Y(1)+T2*Y(2) F(2) = T2*Y(1)+T1*Y(2) RETURN 10 F(1) = 1.+B F(2) = 1.-B RETURN END SUBROUTINE ADBASH (H,X,N,YIP1,YI,YIM1,YIM2,YIM3,F,G) C C--------------------------------------------------------------------- C C Subroutine for performing one step in the solution of a C system of first-order differential equations by the C Adams-Bashford method. C C H is the step size C X is the initial value of the independent variable C N is the number of dependent variables (or equations) C YIP1 is an array of values to be calculated C YI is an array of values from the previous time step C YIM1 is an array of values from two time steps before C YIM2 is an array of values from three time steps before C YIM3 is an array of values from four time steps before C F,G are arrays dimensioned at least N to be supplied by C the calling program. They are scratch arrays used for C intermediate calculations. C C The initial values of the arrays YI, YIM1, YIM2, and YIM3 C must be supplied by the user by some other method. The C subroutine KUTSIM, which requires only the starting value C can be used to compute the starting values at the next C three time steps. C C The arrays YI, YIM1, YIM2, YIM3 are updated one time step C before exit from this subroutine. C C The user must supply a function evaluation routine C SUBROUTINE FUNC (F, X, Y) C that returns the values of the N derivative functions in C the vector F. C C Written by Robert Walraven 9 Oct 81 C Applied Science C Univ. of Calif., Davis 95616 C (916) 752-0360 C C--------------------------------------------------------------------- C DIMENSION YIP1(1),YI(1),YIM1(1),YIM2(1),YIM3(1),F(1),G(1) C H24 = H/24. CALL FUNC (F, X, YI) DO 10 I=1,N G(I) = 55. * F(I) 10 YIP1(I) = 19. * F(I) CALL FUNC (F, X-H, YIM1) DO 20 I=1,N G(I) = G(I) - 59. * F(I) 20 YIP1(I) = YIP1(I) - 5. * F(I) CALL FUNC (F, X-2.*H, YIM2) DO 30 I=1,N G(I) = G(I) + 37. * F(I) 30 YIP1(I) = YIP1(I) + F(I) CALL FUNC (F, X-3.*H, YIM3) DO 40 I=1,N 40 G(I) = YI(I) + H24 * (G(I)-9.*F(I)) CALL FUNC (F, X+H, G) DO 50 I=1,N YIP1(I) = YI(I) + H24 * (YIP1(I)+9.*F(I)) 50 YIP1(I) = (19.*G(I)+251.*YIP1(I))/270. DO 60 I=1,N YIM3(I) = YIM2(I) YIM2(I) = YIM1(I) YIM1(I) = YI(I) 60 YI(I) = YIP1(I) RETURN END SUBROUTINE KUTSIM (H, X, N, Y1, Y2, F, G) C C--------------------------------------------------------------------- C C Subroutine for performing one step in the solution of a C system of first-order differential equations by the C Kutta-Simpson method. C C H is the step size C X is the initial value of the independent variable C N is the number of dependent variables (or equations) C Y1 is an array of initial values of the dependent variables C Y2 is an array of final values of the dependent variables C F,G are arrays dimensioned at least N to be supplied by C the calling program. They are scratch arrays used for C intermediate calculations. C C The user must supply a function evaluation routine C SUBROUTINE FUNC (F, X, Y) C that returns the values of the N derivative functions in C the vector F. C C Written by Robert Walraven 9 Oct 81 C Applied Science C Univ. of Calif., Davis 95616 C (916) 752-0360 C C--------------------------------------------------------------------- C DIMENSION Y1(1),Y2(1),F(1),G(1) C CALL FUNC (F,X,Y1) DO 10 I=1,N G(I) = H * F(I) Y2(I) = G(I) 10 G(I) = Y1(I) + 0.5 * G(I) CALL FUNC (F, X+0.5*H, G) DO 20 I=1,N G(I) = H * F(I) Y2(I) = Y2(I) + 2. * G(I) 20 G(I) = Y1(I) + 0.5 * G(I) CALL FUNC (F, X+0.5*H, G) DO 30 I=1,N G(I) = H * F(I) Y2(I) = Y2(I) + 2. * G(I) 30 G(I) = Y1(I) + G(I) CALL FUNC (F, X+H, G) DO 40 I=1,N 40 Y2(I) = Y1(I) + (Y2(I)+H*F(I))/6. RETURN END