C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- PROGRAM RUNGE C Program: Numerically solves arbitrary ordinary differential equations C using the general coupled 4th-order Runge Kutta problem written as C C u'(1) = f1(x,u(1),u(2),...,u(n)) C u'(2) = f2(x,u(1),u(2),...,u(n)) C . . . C . . . C . . . C u'(n) = fn(x,u(1),u(2),...,u(n)) C C for the dependent variables u(1),u(2),...,u(n). C C The code is written with n=10, meaning that differential equations C up to and including 10th order can be numerically solved with this C program. C C Example (the general 2nd-order differential equation): C C y''(x) + p(x)y'(x) + q(x)y(x) = r(x) (1) C C In order, define C C u(1) = y'(x) C u(2) = y(x) (2) C C (In general, for an Nth-order differential equation y(x) is C always u(N), y'(x) is u(N-1), y''(x) is u(N-2), etc.) C C From (2), the above differential equation (1) becomes C C u'(1) + p(x)u(1) + q(x)u(2) = r(x). C C From (2) it also follows trivially that u'(2) = u(1). C C Then C u'(1) = r(x) - p(x)u(1) - q(x)u(2) C u'(2) = u(1) (3) C C The set of two 1st-order coupled differential equations in (3) can C then be solved together using this general Runge Kutta program after C proper editing of functions f1 and f2 below (you must make sure that C the other function return values are zero - i.e., F3=0, F4=0, ..., C F10=0). The numerical result for u(2) after the calculations is C identically the desired dependent variable y(x). C C This program shows an example of numerically solving a 3rd-order C differential equation: C C y'''(x) + cos(x)y'(x) + xy = sin(x) - x*x (4) C C Using the same ordering criterion as with the above 2nd-order C differential equation, C C u(1) = y''(x) C u(2) = y'(x) C u(3) = y(x), C C it follows that C C u'(1) = sin(x) - x*x - cos(x)u(2) - x*u(3) C u'(2) = u(1) C u'(3) = u(2) (5) C C and the final desired dependent variable y(x) is u(3). This program C prints out a table of incremental values of independent variable x, C u(1), u(2), and u(3) to an ASCII file named "runge.out", and also to C the screen. In the program you must edit the beginning and ending C value of the independent variable x, and the number of incremental C steps the Runge Kutta integration between these beginning and ending C values of x. You must also state the number of values skipped in C the printed table. I.e., with number of steps = 10000, choosing the C number skipped at 1000 means that only 10 lines are printed in the C table. Choosing number skipped = 0 means that 10000 lines (i.e., C all iterative integration steps) are printed to the table. C C With every problem you must also state the initial boundary C conditions. For this example program involving a 3rd-order C differential equation, the beginning boundary conditions were C chosen as follows (A = beginning value for x): C C y(A) = 1 <== u(3) @ x=A C y'(A) = 2 <== u(2) @ x=A C y''(A) = 3 <== u(1) @ x=A C C Note: For 1st, 2nd, and 3rd-order differential equations the printed C ASCII table provides everything you need to know for x value, C derivatives, and the final y(x). y(x) in the printed table is u(1) C for 1st order, u(2) for 2nd order, and u(3) for 3rd order. For 1st C order, ignore u(2) and u(3) in the table, and for 2nd order ignore C u(3). For differential equations with order greater than 3 you will C need to edit the WRITE statements. As example, for a 6th-order C differential equation you might change u(1), u(2), u(3) in the WRITE C statements to u(4), u(5), u(6) to list y''(x), y'(x), and y(x), C respectively. C C STEPS: C 1. Make a backup copy of this original program. C 2. Go through above steps as shown for 2nd and 3rd-order C differential equations to arrive at your set of C 1st-order coupled differential equations for U'(1), C etc. (e.g., like (3) or (5) above). C 3. Edit this program and change functions F1, F2, etc. C accordingly. Make sure that return values for C function numbers greater than the order of your C differential equation are set to zero (e.g., for a C 4th-order differential equation, F5, F6, ..., F10 must C all have zero return value). C 4. Edit this program and accordingly change parameters A, C B, NUM, SKIP, and the initial boundary conditions y(A), C y'(A), etc. C 5. For differential equations with order greater than C three, you will need to edit the WRITE statements in C this program (see discussion in previous paragraph). C C Author: Jerry R. Ziemke PARAMETER(PI=3.141592653589793) INTEGER COUNT,SKIP REAL U(10),K0(10),K1(10),K2(10),K3(10) !*** MUST EDIT THE FOLLOWING (Number of incremental steps (NUM), !initial value of independent variable x (A), ending value of !x (B), and number of iterative steps skipped in the printed !table (SKIP)): NUM=10000 ! <== Number of incremental steps between A and B A=-5 ! <== Beginning value for independent variable X B=5 ! <== Ending value for independent variable X SKIP=10 ! <== Print values only every "SKIP" # of times !*** MUST EDIT THE FOLLOWING (initial condition values): U(1)=1.0 ! <== Initial value y''(A) for 3rd-order example U(2)=2.0 ! <== Initial value y'(A) for 3rd-order example U(3)=3.0 ! <== Initial value y(A) for 3rd-order example OPEN(12,FORM='FORMATTED',STATUS='UNKNOWN',FILE='runge.out') WRITE(12,*) & ' X Y''''(X) Y''(X) Y(X)' WRITE(12,*) & '-------------------------------------------------------' WRITE(*,*) & ' X Y''''(X) Y''(X) Y(X)' WRITE(*,*) & '-------------------------------------------------------' H=(B-A)/NUM COUNT=0 DO I=1,NUM+1 X=A+(I-1)*H COUNT=COUNT+1 IF ((COUNT.EQ.(SKIP+1)).OR.(I.EQ.1).OR.(SKIP.EQ.0)) THEN COUNT=0 WRITE(*,'(4E14.5)') X, U(1), U(2), U(3) WRITE(12,'(4E14.5)') X, U(1), U(2), U(3) ENDIF U1=U(1) U2=U(2) U3=U(3) U4=U(4) U5=U(5) U6=U(6) U7=U(7) U8=U(8) U9=U(9) U10=U(10) K0(1)=H*F1(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K0(2)=H*F2(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K0(3)=H*F3(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K0(4)=H*F4(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K0(5)=H*F5(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K0(6)=H*F6(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K0(7)=H*F7(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K0(8)=H*F8(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K0(9)=H*F9(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K0(10)=H*F10(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) U1=U1+K0(1)*0.5 U2=U2+K0(2)*0.5 U3=U3+K0(3)*0.5 U4=U4+K0(4)*0.5 U5=U5+K0(5)*0.5 U6=U6+K0(6)*0.5 U7=U7+K0(7)*0.5 U8=U8+K0(8)*0.5 U9=U9+K0(9)*0.5 U10=U10+K0(10)*0.5 K1(1)=H*F1(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K1(2)=H*F2(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K1(3)=H*F3(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K1(4)=H*F4(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K1(5)=H*F5(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K1(6)=H*F6(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K1(7)=H*F7(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K1(8)=H*F8(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K1(9)=H*F9(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K1(10)=H*F10(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) U1=U1+(K1(1)-K0(1))*0.5 U2=U2+(K1(2)-K0(2))*0.5 U3=U3+(K1(3)-K0(3))*0.5 U4=U4+(K1(4)-K0(4))*0.5 U5=U5+(K1(5)-K0(5))*0.5 U6=U6+(K1(6)-K0(6))*0.5 U7=U7+(K1(7)-K0(7))*0.5 U8=U8+(K1(8)-K0(8))*0.5 U9=U9+(K1(9)-K0(9))*0.5 U10=U10+(K1(10)-K0(10))*0.5 K2(1)=H*F1(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K2(2)=H*F2(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K2(3)=H*F3(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K2(4)=H*F4(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K2(5)=H*F5(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K2(6)=H*F6(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K2(7)=H*F7(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K2(8)=H*F8(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K2(9)=H*F9(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K2(10)=H*F10(X+H*0.5,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) U1=U1+K2(1)-K1(1)*0.5 U2=U2+K2(2)-K1(2)*0.5 U3=U3+K2(3)-K1(3)*0.5 U4=U4+K2(4)-K1(4)*0.5 U5=U5+K2(5)-K1(5)*0.5 U6=U6+K2(6)-K1(6)*0.5 U7=U7+K2(7)-K1(7)*0.5 U8=U8+K2(8)-K1(8)*0.5 U9=U9+K2(9)-K1(9)*0.5 U10=U10+K2(10)-K1(10)*0.5 K3(1)=H*F1(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K3(2)=H*F2(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K3(3)=H*F3(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K3(4)=H*F4(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K3(5)=H*F5(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K3(6)=H*F6(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K3(7)=H*F7(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K3(8)=H*F8(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K3(9)=H*F9(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) K3(10)=H*F10(X+H,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) U(1)=U(1)+(K0(1)+2*K1(1)+2*K2(1)+K3(1))/6 U(2)=U(2)+(K0(2)+2*K1(2)+2*K2(2)+K3(2))/6 U(3)=U(3)+(K0(3)+2*K1(3)+2*K2(3)+K3(3))/6 U(4)=U(4)+(K0(4)+2*K1(4)+2*K2(4)+K3(4))/6 U(5)=U(5)+(K0(5)+2*K1(5)+2*K2(5)+K3(5))/6 U(6)=U(6)+(K0(6)+2*K1(6)+2*K2(6)+K3(6))/6 U(7)=U(7)+(K0(7)+2*K1(7)+2*K2(7)+K3(7))/6 U(8)=U(8)+(K0(8)+2*K1(8)+2*K2(8)+K3(8))/6 U(9)=U(9)+(K0(9)+2*K1(9)+2*K2(9)+K3(9))/6 U(10)=U(10)+(K0(10)+2*K1(10)+2*K2(10)+K3(10))/6 ENDDO CLOSE(12) STOP END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F1(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F1=SIN(X)-X*X+U2*COS(X)-X*U3 RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F2(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F2=U1 RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F3(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F3=U2 RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F4(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F4=0 RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F5(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F5=0 RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F6(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F6=0 RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F7(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F7=0 RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F8(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F8=0 RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F9(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F9=0 RETURN END C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-- FUNCTION F10(X,U1,U2,U3,U4,U5,U6,U7,U8,U9,U10) F10=0 RETURN END