'******************************************************************** '* Differential equations with p variables of order 1 * '* by Runge-Kutta method of order 4 * '* ---------------------------------------------------------------- * '* Reference: "Analyse en Turbo Pascal versions 5.5 et 6.0 by Marc * '* DUCAMP et Alain REVERCHON - Eyrolles, Paris 1991" * '* * '* Basic version 1.1 By J-P Moreau, Paris * '* (www.jpmoreau.fr) * '* ---------------------------------------------------------------- * '* SAMPLE RUN: * '* * '* Example #1: integrate system of equations from x=0 to x=3: * '* y1' = y2 + y3 - 3*y1 * '* y2' = y1 + y3 - 3*y2 * '* y3' = y1 + y2 - 3*y3 * '* with initial conditions: y1(0)=1, y2(0)=2 and y3(0)=-1 * '* * '* * '* DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * '* of type yi' = f(y1,y2,...,yn), i=1..n * '* * '* number of variables: 3 * '* begin value x : 0 * '* end value x : 3 * '* y1 value at x0 : 1 * '* y2 value at x0 : 2 * '* y3 value at x0 : -1 * '* number of points : 7 * '* finesse : 30 * '* * '* X Y1 Y2 Y3 * '* -------------------------------------------------------- * '* 0.000000 1.000000 2.000000 -1.000000 * '* 0.500000 0.449466 0.584801 0.178795 * '* 1.000000 0.251358 0.269674 0.214727 * '* 1.500000 0.149580 0.152058 0.144622 * '* 2.000000 0.090335 0.090671 0.089664 * '* 2.500000 0.054738 0.054784 0.054648 * '* 3.000000 0.033193 0.033200 0.033181 * '* -------------------------------------------------------- * '* * '* Example #2: integrate system of equations from x=0 to PI*Sqrt(2) * '* y1' = y2 * '* y2' = -4*y1 - 3*y3 * '* y3' = y4 * '* y4' = -8*y1 - 2*y3 * '* with initial conditions: y1(0)=3, y2(0)=0, y3(0)=4, y4(0)=0 * '* * '* * '* DIFFERENTIAL EQUATION WITH P VARIABLE OF ORDER 1 * '* of type yi' = f(y1,y2,...,yn), i=1..n * '* * '* number of variables: 4 * '* begin value x : 0 * '* end value x : 4.442883 * '* y1 value at x0 : 3 * '* y2 value at x0 : 0 * '* y3 value at x0 : 4 * '* y4 value at x0 : 0 * '* number of points : 9 * '* finesse : 30 * '* * '* X Y1 Y2 Y3 Y4 * '* -------------------------------------------------------- * '* 0.000000 3.000000 0.000000 4.000000 0.000000 * '* 0.555360 0.000000 -8.485281 0.000000 -11.313708 * '* 1.110721 -3.000000 -0.000001 -4.000000 -0.000002 * '* 1.666081 -0.000001 8.485281 -0.000001 11.313708 * '* 2.221442 3.000000 0.000003 4.000000 0.000003 * '* 2.776802 0.000001 -8.485281 0.000002 -11.313708 * '* 3.332162 -3.000000 -0.000004 -4.000000 -0.000005 * '* 3.887523 -0.000002 8.485281 -0.000002 11.313708 * '* 4.442883 3.000000 0.000005 4.000000 0.000007 * '* -------------------------------------------------------- * '* * '* Release 1.1: Added example #2 (Sept. 2008). * '******************************************************************** DefInt I-N DefDbl A-H, O-Z 'ifi,ip : INTEGER Dim y(10), yi(10), t1(50), t2(50), t3(50), t4(50) Cls Print Print " DIFFERENTIAL EQUATIONS WITH P VARIABLES OF ORDER 1" Print " of type yi' = f(y1,y2,...,yn), i=1..n" Print PRINT " number of variables : "; : INPUT ip Print PRINT " begin value x : "; : INPUT xi PRINT " end value x : "; : INPUT xf For i = 0 To ip - 1 PRINT " y"; i + 1; " value at x0 : "; : INPUT yi(i) Next i PRINT " number of points : "; : INPUT m PRINT " finesse : "; : INPUT ifi 'call subroutine eqdifp GoSub 2000 End 'Example #1: y1'=y2+y3-3y1, y2'=y1+y3-3y2, y3'=y1+y2-3y3 '1000 'FUNCTION fp ' IF k = 0 THEN ' fp = y(1) + y(2) - 3# * y(0) ' ELSEIF k = 1 THEN ' fp = y(0) + y(2) - 3# * y(1) ' ELSEIF k = 2 THEN ' fp = y(0) + y(1) - 3# * y(2) ' ELSE ' fp = 0# ' END IF 'RETURN 'Example #2: y1'=y2, y2'=-4y1-3y3, y'3=y4, y'4=-8y1-2y3 1000 'FUNCTION fp If k = 0 Then fp = y(1) ElseIf k = 1 Then fp = -4# * y(0) - 3# * y(2) ElseIf k = 2 Then fp = y(3) ElseIf k = 3 Then fp = -8# * y(0) - 2# * y(2) Else fp = 0# 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) t3(1) = yi(2) t4(1) = yi(3) 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) t3(i + 1) = z(2) t4(i + 1) = z(3) Next i 'Display(t1,t2,t3,t4,m,p,xi,xf) GoSub 3000 Return 3000 'subroutine Display h = (xf - xi) / (m - 1) X = xi - h Print Print " X"; For i = 1 To ip + 1 Print " Y"; i; Next i Print Print "--------------------------------------------------------" For i = 1 To m X = X + h If ip = 1 Then Print USING; " ##.###### ##.###### ##.######"; X; t1(i); t2(i) ElseIf ip = 2 Then Print USING; " ##.###### ##.###### ##.###### ##.######"; X; t1(i); t2(i); t3(i) ElseIf ip > 2 Then Print USING; " ##.###### ##.###### ##.###### ##.###### ###.######"; X; t1(i); t2(i); t3(i); t4(i) End If Next i Print "--------------------------------------------------------" Return 'End of file teqdifp.bas
Wednesday, September 20, 2017
Runge-Kutta method of order 4
Subscribe to:
Post Comments (Atom)
No comments:
Post a Comment