'******************************************************************** '* Differential equations with p variables of order 1 * '* by Runge-Kutta method of order 4 * '* (Case of a Mass Pendulum) * '* ---------------------------------------------------------------- * '* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 By Marc * '* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" * '* [BIBLI 03]. * '* * '* Basic Release By J-P Moreau, Paris. * '* (www.jpmoreau.fr) * '* ---------------------------------------------------------------- * '* EXPLANATION: * '* Let us consider the following Mass Pendulum: * '* * '* . O * '* /| * '* d /_|theta * '* / | * '* / | * '* | * '* * '* The deviation angle theta is determined by the differential * '* equation: * '* theta" = -K sin(theta) - K' theta * '* where: * '* K is a coefficient depending on the pendulum * '* characteristics (K=Mgd/J), * '* K' is another coefficient for fluid resistance. * '* * '* Moreover, one must know the initial conditions theta(0) and * '* theta'(0). * '* * '* This differential equation will be replaced by the system with * '* two variables: * '* y1' = y2 * '* y'2 = -K*sin(y1) -K'*y2 * '* * '* (this allows calculating theta and theta' at each time step). * '* We use a Runge-Kutta method of 4th order to solve the system * '* with respect to time. Here, K=(pi^2/4) and K'=0.1 (see unit * '* Eqdifp.pas). * '* ---------------------------------------------------------------- * '* SAMPLE RUNS: * '* * '* Example1: integrate system of equations from t=0 to t=12: * '* y1' = y2 * '* y2' = -(pi^2/4)*sin(y1) - 0.1*y2 * '* with initial conditions: y1(0)=0.25, y2(0)=0 * '* * '* * '* DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * '* of type yi' = f(y1,y2,...,yn), i=1..n * '* (Case of a Mass Pendulum) * '* * '* Begin value time : 0 * '* End value time : 12 * '* theta value at t0 : 0.25 * '* theta' value at t0 : 0 * '* Number of points : 25 * '* Finesse : 10 * '* * '* Time Theta Theta' * '* -------------------------------------- * '* 0.000000 0.250000 0.000000 * '* 0.500000 0.178615 -0.268789 * '* 1.000000 0.009250 -0.372792 * '* 1.500000 -0.157126 -0.259240 * '* 2.000000 -0.226032 -0.004518 * '* 2.500000 -0.163373 0.240320 * '* 3.000000 0.010873 0.337319 * '* 3.500000 0.140388 0.237242 * '* 4.000000 0.204362 0.007531 * '* 4.500000 0.149147 -0.215220 * '* 5.000000 0.011745 -0.305188 * '* 5.500000 -0.125654 -0.216683 * '* 6.000000 -0.184778 -0.009451 * '* 6.500000 -0.135956 0.193012 * '* 7.000000 -0.012089 0.276100 * '* 7.500000 0.112634 0.197596 * '* 8.000000 0.167078 0.010581 * '* 8.500000 0.123784 -0.173302 * '* 9.000000 0.012066 -0.249776 * '* 9.500000 -0.101090 -0.179966 * '* 10.000000 -0.151083 -0.011147 * '* 10.500000 -0.112595 0.155759 * '* 11.000000 -0.011793 0.225958 * '* 11.500000 0.090824 0.163747 * '* 12.000000 0.136627 0.011313 * '* * '* Example2: integrate system of equations from t=0 to t=12: * '* y1' = y2 * '* y2' = -(pi^2/4)*sin(y1) - 0.1*y2 * '* with initial conditions: y1(0)=3.14, y2(0)=0 * '* * '* * '* DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * '* of type yi' = f(y1,y2,...,yn), i=1..n * '* (Case of a Mass Pendulum) * '* * '* Begin value time : 0 * '* End value time : 40 * '* theta value at t0 : 3.14 * '* theta' value at t0 : 0 * '* Number of points : 41 * '* Finesse : 10 * '* * '* Time Theta Theta' * '* -------------------------------------- * '* 0.000000 3.140000 0.000000 * '* 1.000000 3.137678 -0.005478 * '* 2.000000 3.124331 -0.026170 * '* 3.000000 3.062686 -0.120013 * '* 4.000000 2.781242 -0.545197 * '* 5.000000 1.576348 -2.134046 * '* 6.000000 -1.105676 -2.312060 * '* 7.000000 -2.175457 0.053644 * '* 8.000000 -1.057882 2.213567 * '* 9.000000 1.294094 1.644802 * '* 10.000000 1.684059 -0.817529 * '* 11.000000 -0.130730 -2.326526 * '* 12.000000 -1.563901 -0.249465 * '* 13.000000 -0.636818 1.933854 * '* 14.000000 1.166974 1.019504 * '* 15.000000 0.988226 -1.309951 * '* 16.000000 -0.727610 -1.446374 * '* 17.000000 -1.088264 0.765154 * '* 18.000000 0.357682 1.576103 * '* 19.000000 1.058234 -0.342549 * '* 20.000000 -0.086719 -1.524313 * '* 21.000000 -0.968792 0.032876 * '* 22.000000 -0.095485 1.389859 * '* 23.000000 0.858794 0.182214 * '* 24.000000 0.210044 -1.230693 * '* 25.000000 -0.748336 -0.322701 * '* 26.000000 -0.276785 1.074664 * '* 27.000000 0.646728 0.407293 * '* 28.000000 0.311067 -0.933079 * '* 29.000000 -0.557369 -0.451786 * '* 30.000000 -0.323904 0.809220 * '* 31.000000 0.480640 0.468557 * '* 32.000000 0.322919 -0.702785 * '* 33.000000 -0.415543 -0.466803 * '* 34.000000 -0.313287 0.612046 * '* 35.000000 0.360579 0.453126 * '* 36.000000 0.298489 -0.534862 * '* 37.000000 -0.314185 -0.432143 * '* 38.000000 -0.280842 0.469135 * '* 39.000000 0.274922 0.407030 * '* 40.000000 0.261876 -0.412985 * '* * '* You may try other examples, such as: * '* theta0=4, theta'0=0 from t=0 to t=20 * '* theta0=0, theta'0=20 from t=0 to t=10 * '* theta0=0, theta'0=5 from t=0 to t=10 * '* theta0=0, theta'0=3.15 from t=0 to t=20 * '* * '******************************************************************** DEFINT I-N DEFDBL A-H, O-Z OPTION BASE 0 'ifi,ip : INTEGER DIM y(10), yi(10), t1(50), t2(50) CLS PRINT PRINT " DIFFERENTIAL EQUATIONS WITH P VARIABLES OF ORDER 1" PRINT " of type yi' = f(y1,y2,...,yn), i=1..n" PRINT " (case of a Mass Pendulum)" PRINT ip = 2 INPUT " begin value time : ", xi INPUT " end value time : ", xf INPUT " Theta at t0 : ", yi(0) INPUT " Theta' at t0 : ", yi(1) INPUT " number of points : ", m INPUT " finesse : ", ifi 'call subroutine eqdifp GOSUB 2000 END 'Example : y1'=y2, y2'=(pi*pi/4)*sin(y1)-0.1*y2 1000 'FUNCTION fp PI = 4# * ATN(1#) IF k = 0 THEN fp = y(1) ELSEIF k = 1 THEN fp = -(PI * PI / 4#) * SIN(y(0)) - .1# * y(1) END IF RETURN '**************************************************************************** '* SOLVING DIFFERENTIAL SYSTEMS WITH P VARIABLES OF ORDER 1 * '* of type yi' = f(y1,y2,...,yn), i=1..n * '* ------------------------------------------------------------------------ * '* INPUTS: * '* m number of points to calculate * '* xi, xf begin, end values of variable x * '* yi table of begin values of functions at xi * '* ip number of independant variables * '* ifi finesse (number of intermediary points) * '* * '* OUTPUTS: * '* t1,t2 real vectors storing the results for first two functions, * '* y1 and y2. * '* ------------------------------------------------------------------------ * '* EXAMPLE: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3 * '* Exact solution : y1 = 1/3 (exp(-4x) + 2 exp(-x)) * '* y2 = 1/3 (4exp(-4x) + 2 exp(-x)) * '* y3 = 1/3 (-5exp(-4x)+ 2 exp(-x)) * '**************************************************************************** 2000 'subroutine Eqdifp DIM ta(10), tb(10), tc(10), td(10), z(10) h = (xf - xi) / ifi / (m - 1) ip = ip - 1 t1(1) = yi(0) t2(1) = yi(1) FOR k = 0 TO ip y(k) = yi(k): z(k) = yi(k) NEXT k FOR i = 1 TO m ni = (i - 1) * ifi - 1 FOR j = 1 TO ifi x = xi + h * (ni + j) FOR k = 0 TO ip y(k) = z(k) NEXT k FOR k = 0 TO ip GOSUB 1000 ta(k) = h * fp NEXT k FOR k = 0 TO ip y(k) = z(k) + ta(k) / 2# NEXT k x = x + h / 2# FOR k = 0 TO ip GOSUB 1000 tb(k) = h * fp NEXT k FOR k = 0 TO ip y(k) = z(k) + tb(k) / 2# NEXT k FOR k = 0 TO ip GOSUB 1000 tc(k) = h * fp NEXT k FOR k = 0 TO ip y(k) = z(k) + tc(k) NEXT k x = x + h / 2# FOR k = 0 TO ip GOSUB 1000 td(k) = h * fp NEXT k FOR k = 0 TO ip z(k) = z(k) + (ta(k) + 2# * tb(k) + 2# * tc(k) + td(k)) / 6# NEXT k NEXT j t1(i + 1) = z(0) t2(i + 1) = z(1) NEXT i 'Sisplay(t1,t2,m,xi,xf) GOSUB 3000 RETURN 3000 'subroutine Display h = (xf - xi) / (m - 1) x = xi - h CLS PRINT PRINT " Time Theta Theta' " PRINT "----------------------------------" FOR i = 1 TO m x = x + h PRINT USING " ##.###### ##.###### ##.######"; x; t1(i); t2(i) NEXT i PRINT RETURN 'End of file pendulum.bas