{$A+,B-,D+,E+,F+,G-,I+,L+,N+,O-,P-,Q-,R-,S+,T+,V+,X+,Y+} unit RKF; { Runge-Kutta-Fehlberg Criado por Rog‚rio Luiz em MODULA-3 Traduzido (e adaptado) para TURBO PASCAL por Forster } INTERFACE TYPE real = double; Vector = ARRAY [0..999] OF REAL; {tipo vetor generico} PVector = ^Vector; {apontador para qualquer vetor} TDerivative = PROCEDURE (t: REAL; n: integer; s, v: PVector); { Dado um tempo "t" e o estado correspondente "s[0..n-1]", calcula a derivada do estado "v[0..n-1]". } TDetectEvent = FUNCTION (t0, t1: REAL; n: integer; s0, s1: PVector): REAL; { Dados dois tempos "t0" e "t1", e os estados correspondentes "s0[0..n-1]" e "s1[1..n-1]", determina o instante "t" no intervalo "[t0 __ t1]" onde ocorre o proximo evento. } TUseState = PROCEDURE (t: REAL; n: integer; s: PVector); { Usa o estado "s[0..n-1]" correspondente ao tempo "t". } P = ^T; T = OBJECT n: integer; {tamanho} constructor init(An: integer); function integrate( t0, t1, dtmin, dtmax, tol: real; s: PVector; derivative: TDerivative; detectEvent: TDetectEvent; useState: TUseState ): real; { Integra a equacao "ds/dt = F(t, s)" do intante "t0" e estado inicial "s" ate o instante "t1", ou ate o proximo evento (o que vier antes). O estado final e' devolvido no proprio "s". Usa uma ou mais iteracoes do metodo RKF, ajustando o passo entre "dtmin" e "dtmax", de modo a conseguir erro estimado menor que "tol" se possivel. Devolve como resultado o tempo em que a integracao parou. } destructor finish; PRIVATE s0, ds0, ds2, ds3, ds4, ds5, ds6: PVector; END; function NEWVECTOR(n: integer): PVector; procedure GETRIDOFVECTOR (n:integer; v:PVector); PROCEDURE CopyState(n: integer; sfr, sto: PVector); IMPLEMENTATION function NEWVECTOR (n: integer): PVector; var v: PVector; begin Getmem(v, n*sizeof(real)); NEWVECTOR := v; end; procedure GETRIDOFVECTOR (n:integer; v:PVector); begin Freemem(v, n*sizeof(real)); end; CONST C12 = 1.0/4.0; Ct2 = C12; C13 = 3.0/32.0; C23 = 9.0/32.0; Ct3 = C13 + C23; C14 = 1932.0/2197.0; C24 = -7200.0/2197.0; C34 = 7296.0 /2197.0 ; Ct4 = C14 + C24 + C34; C15 = 439.0 /216.0 ; C25 = -8.0 ; C35 = 3680.0 /513.0 ; C45 = -845.0 /4104.0 ; Ct5 = C15 + C25 + C35 + C45; C16 = -8.0 /27.0 ; C26 = 2.0 ; C36 = -3544.0 /2565.0 ; C46 = 1859.0 /4104.0 ; C56 = -11.0 /40.0 ; Ct6 = C16 + C26 + C36 + C46 + C56; Ce1 = 1.0 /360.0 ; Ce3 = -128.0 /4275.0 ; Ce4 = -2197.0 /75240.0 ; Ce5 = 1.0 /50.0 ; Ce6 = 2.0 /55.0 ; Cv1 = 25.0 /216.0 ; Cv3 = 1408.0 /2565.0 ; Cv4 = 2197.0 /4104.0 ; Cv5 = -1.0 /5.0 ; Constructor T.Init(An: integer); BEGIN n := An; s0 := NEWVector(n); ds0 := NEWVector(n); ds2 := NEWVector(n); ds3 := NEWVector(n); ds4 := NEWVector(n); ds5 := NEWVector(n); ds6 := NEWVector(n); END {Init}; function Step( t0, t1: real; n: integer; s0, ds0: PVector; ds2, ds3, ds4, ds5, ds6, s1: PVector; { Areas de trabalho } derivative: TDerivative ): real; { Calcula uma aproximacao de quarta ordem para a trajetoria do sistema entre "t0" e "t1", dado o estado inicial "s0" e a derivada correspondente "ds0". Usa o metodo RK com um unico passo. devolve o erro estimado. Para obter o estado em "t1", chame "FinalState(t0, t1, n, s0, ds0, ds3, ds4, ds5, s1)". } var i: integer; t, tf, dt, rds, ards, error: real; BEGIN t:=t0; dt := t1 - t0; FOR i := 0 TO n-1 DO begin s1^[i] := s0^[i] + dt*C12*ds0^[i]; END; derivative(t + Ct2*dt, n, s1, ds2); FOR i := 0 TO n-1 DO begin s1^[i] := s0^[i] + dt*(C13*ds0^[i] + C23*ds2^[i]); END; derivative(t + Ct3*dt, n, s1, ds3); FOR i := 0 TO n-1 DO begin s1^[i] := s0^[i] + dt*(C14*ds0^[i] + C24*ds2^[i] + C34*ds3^[i]); END; derivative(t + Ct4*dt, n, s1, ds4); FOR i := 0 TO n-1 DO begin s1^[i] := s0^[i] + dt*( C15*ds0^[i] + C25*ds2^[i] + C35*ds3^[i] + C45*ds4^[i] ); END; derivative(t + Ct5*dt, n-1, s1, ds5); FOR i := 0 TO n-1 DO begin s1^[i] := s0^[i] + dt*( C16*ds0^[i] + C26*ds2^[i] + C36*ds3^[i] + C46*ds4^[i] + C56*ds5^[i] ); END; derivative(t + Ct6*dt, n, s1, ds6); error := 0.0; FOR i := 0 TO n-1 DO begin rds := Ce1*ds0^[i] + Ce3*ds3^[i] + Ce4*ds4^[i] + Ce5*ds5^[i] + Ce6*ds6^[i]; ards := ABS(rds); IF ards > error THEN error := ards ; END; Step := error; END {Step}; PROCEDURE FinalState( t0, t1: real; n: integer; s0, ds0, ds3, ds4, ds5, s1: PVector ); { Calcula o estado final "s1" da integracao executada por "Step", dados o estado inicial "s0", sua derivada "ds0", e as derivadas auxiliares "ds3, ds4, ds5" calculadas por "Step". } var i:integer; BEGIN FOR i := 0 TO n-1 DO begin s1^[i] := s0^[i] + (t1-t0)*(Cv1*ds0^[i] + Cv3*ds3^[i] + Cv4*ds4^[i] + Cv5*ds5^[i]); END END {FinalState}; PROCEDURE CopyState(n: integer; sfr, sto: PVector); BEGIN MOVE(sfr^, sto^, n*sizeof(real)); END {CopyState}; PROCEDURE Updatedt(var dt: real; error, tol, dtmin, dtmax: real); var delta:real; BEGIN IF error > 0.0 THEN begin delta := 0.840896415 * sqrt(sqrt(tol/error)); IF delta <= 0.1 THEN dt := 0.1 * dt ELSE IF delta >= 4.0 THEN dt := 4.0 * dt ELSE dt := delta * dt; END ELSE dt := 4.0 * dt; IF dt > dtmax THEN dt := dtmax ELSE IF dt < dtmin THEN dt := dtmin; END {Updatedt}; function T.Integrate( t0, t1, dtmin, dtmax, tol: real; s: PVector; derivative: TDerivative; detectEvent: TDetectEvent; useState: TUseState ): real; VAR ok: BOOLEAN; t, tf, dt, tev, error: REAL; BEGIN t := t0; dt := sqrt(dtmax*dtmin); while t < t1 do begin CopyState(n, s, s0); useState(t, n, s0); derivative(t, n, s0, ds0); { Tenta avancar de "t" para algum instante posterior "tf": } repeat ok := false; tf := t + dt; IF tf > t1 then begin tf := t1; dt := t1 - t end; error := Step(t, tf, n, s0, ds0, ds2, ds3, ds4, ds5, ds6, s, derivative); IF (error <= tol) OR (dt <= dtmin) THEN BEGIN FinalState(t, tf, n, s0, ds0, ds3, ds4, ds5, s); tev := detectEvent(t, tf, n, s0, s); IF (tev > t) and (tev < tf) THEN t1 := tev ELSE ok := true END; UpdateDT(dt, error, tol, dtmin, dtmax) until ok; t := tf; end; Integrate := t1 END {Integrate}; destructor T.finish; begin GETRIDOFVECTOR(n,s0); GETRIDOFVECTOR(n,ds0); GETRIDOFVECTOR(n,ds2); GETRIDOFVECTOR(n,ds3); GETRIDOFVECTOR(n,ds4); GETRIDOFVECTOR(n,ds5); GETRIDOFVECTOR(n,ds6); end; BEGIN END {RKF}.