{ This file is part of the Numlib package. Copyright (c) 1986-2000 by Kees van Ginneken, Wil Kortsmit and Loek van Reij of the Computational centre of the Eindhoven University of Technology FPC port Code by Marco van de Voort (marco@freepascal.org) documentation by Michael van Canneyt (Michael@freepascal.org) Solve first order starting value differential eqs, and sets of first order starting value differential eqs, Both versions are not suited for stiff differential equations See the file COPYING.FPC, included in this distribution, for details about the copyright. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. **********************************************************************} Unit ode; {$I DIRECT.INC} interface uses typ; {Solve first order, starting value, differential eqs, Calc y(b) for dy/dx=f(x,y) and y(a)=ae} Procedure odeiv1(f: rfunc2r; a, ya: ArbFloat; Var b, yb: ArbFloat; ae: ArbFloat; Var term: ArbInt); { The same as above, for a set of equations. ya and yb are vectors} Procedure odeiv2(f: oderk1n; a: ArbFloat; Var ya, b, yb: ArbFloat; n: ArbInt; ae: ArbFloat; Var term: ArbInt); implementation Procedure odeiv1(f: rfunc2r; a, ya: ArbFloat; Var b, yb: ArbFloat; ae: ArbFloat; Var term: ArbInt); Var last, first, reject, goon : boolean; x, y, d, h, xl, yl, int, hmin, absh,k0, k1, k2, k3, k4, k5, discr, tol, mu, mu1, fh, hl : ArbFloat; Begin x := a; y := ya; d := b-a; yb := y; term := 1; If ae <= 0 Then Begin term := 3; exit End; If d <> 0 Then Begin xl := x; yl := y; h := d/4; absh := abs(h); int := abs(d); hmin := int*1e-6; ae := ae/int; first := true; goon := true; while goon Do Begin absh := abs(h); If absh < hmin Then Begin If h>0 Then h := hmin Else h := -hmin; absh := hmin End; If (h >= b-xl) = (h >= 0) Then Begin last := true; h := b-xl; absh := abs(h) End Else last := false; x := xl; y := yl; k0 := f(x,y)*h; x := xl+h*2/9; y := yl+k0*2/9; k1 := f(x,y)*h; x := xl+h/3; y := yl+(k0+k1*3)/12; k2 := f(x,y)*h; x := xl+h/2; y := yl+(k0+k2*3)/8; k3 := f(x,y)*h; x := xl+h*0.8; y := yl+(k0*53-k1*135+k2*126+k3*56)/125; k4 := f(x,y)*h; If last Then x := b Else x := xl+h; y := yl+(k0*133-k1*378+k2*276+k3*112+k4*25)/168; k5 := f(x,y)*h; discr := abs(21*(k0-k3)-162*(k2-k3)-125*(k4-k3)+42*(k5-k3))/14; tol := absh*ae; mu := 1/(1+discr/tol)+0.45; reject := discr > tol; If reject Then Begin If absh <= hmin Then Begin b := xl; yb := yl; term := 2; exit End; h := mu*h End Else Begin If first Then Begin first := false; hl := h; h := mu*h End Else Begin fh := mu*h/hl+mu-mu1; hl := h; h := fh*h End; mu1 := mu; y := yl+(-k0*63+k1*189-k2*36-k3*112+k4*50)/28; k5 := f(x,y)*hl; y := yl+(k0*35+k2*162+k4*125+k5*14)/336; If b <> x Then Begin xl := x; yl := y End Else Begin yb := y; goon := false End End {not reject} End; {while} End {d<>0} End; {odeiv1} Procedure odeiv2(f: oderk1n; a: ArbFloat; Var ya, b, yb: ArbFloat; n: ArbInt; ae: ArbFloat; Var term: ArbInt); Var pya, pyb, yl, k0, k1, k2, k3, k4, k5, y : ^arfloat1; i, jj, ns : ArbInt; last, first, reject, goon : boolean; x, xl, hmin, int, hl, absh, fhm, discr, tol, mu, mu1, fh, d, h : ArbFloat; Begin If (ae <= 0) Or (n < 1) Then Begin term := 3; exit End; ns := n*sizeof(ArbFloat); pya := @ya; pyb := @yb; move(pya^[1], pyb^[1], ns); term := 1; getmem(yl, ns); getmem(k0, ns); getmem(k1, ns); getmem(k2, ns); getmem(k3, ns); getmem(k4, ns); getmem(k5, ns); getmem(y, ns); x := a; d := b-a; move(pya^[1], y^[1], ns); If d <> 0 Then Begin xl := x; move(y^[1], yl^[1], ns); h := d/4; absh := abs(h); int := abs(d); hmin := int*1e-6; hl := ae; ae := ae/int; first := true; goon := true; while goon Do Begin absh := abs(h); If absh < hmin Then Begin If h > 0 Then h := hmin Else h := -hmin; absh := hmin End; If (h >= b-xl) = (h >= 0) Then Begin last := true; h := b-xl; absh := abs(h) End Else last := false; x := xl; move(yl^[1], y^[1], ns); f(x, y^[1], k0^[1]); For i:=1 To n Do k0^[i] := k0^[i]*h; x := xl+h*2/9; For jj:=1 To n Do y^[jj] := yl^[jj]+k0^[jj]*2/9; f(x, y^[1], k1^[1]); For i:=1 To n Do k1^[i] := k1^[i]*h; x := xl+h/3; For jj:=1 To n Do y^[jj] := yl^[jj]+(k0^[jj]+k1^[jj]*3)/12; f(x, y^[1], k2^[1]); For i:=1 To n Do k2^[i] := k2^[i]*h; x := xl+h/2; For jj:=1 To n Do y^[jj] := yl^[jj]+(k0^[jj]+k2^[jj]*3)/8; f(x, y^[1], k3^[1]); For i:=1 To n Do k3^[i] := k3^[i]*h; x := xl+h*0.8; For jj:=1 To n Do y^[jj] := yl^[jj]+ (k0^[jj]*53-k1^[jj]*135+k2^[jj]*126+k3^[jj]*56)/125; f(x, y^[1], k4^[1]); For i:=1 To n Do k4^[i] := k4^[i]*h; If last Then x := b Else x := xl+h; For jj:=1 To n Do y^[jj] := yl^[jj]+(k0^[jj]*133-k1^[jj]*378+k2^[jj]*276+ k3^[jj]*112+k4^[jj]*25)/168; f(x, y^[1], k5^[1]); For i:=1 To n Do k5^[i] := k5^[i]*h; reject := false; fhm := 0; tol := absh*ae; For jj:=1 To n Do Begin discr := abs((k0^[jj]-k3^[jj])*21-(k2^[jj]-k3^[jj])*162- (k4^[jj]-k3^[jj])*125+(k5^[jj]-k3^[jj])*42)/14; reject := (discr > tol) Or reject; fh := discr/tol; If fh > fhm Then fhm := fh End; {jj} mu := 1/(1+fhm)+0.45; If reject Then Begin If absh <= hmin Then Begin b := xl; move(yl^[1], pyb^[1], ns); term := 2; freemem(yl, ns); freemem(k0, ns); freemem(k1, ns); freemem(k2, ns); freemem(k3, ns); freemem(k4, ns); freemem(k5, ns); freemem(y, ns); exit End; h := mu*h End Else Begin If first Then Begin first := false; hl := h; h := mu*h End Else Begin fh := mu*h/hl+mu-mu1; hl := h; h := fh*h End; mu1 := mu; For jj:=1 To n Do y^[jj] := yl^[jj]+(-k0^[jj]*63+k1^[jj]*189 -k2^[jj]*36-k3^[jj]*112+k4^[jj]*50)/28; f(x, y^[1], k5^[1]); For i:=1 To n Do k5^[i] := k5^[i]*hl; For jj:=1 To n Do y^[jj] := yl^[jj]+(k0^[jj]*35+k2^[jj]*162+k4^[jj]*125 +k5^[jj]*14)/336; If b <> x Then Begin xl := x; move(y^[1], yl^[1], ns) End Else Begin move(y^[1], pyb^[1], ns); goon := false End End {not reject} End {while} End; {d<>0} freemem(yl, ns); freemem(k0, ns); freemem(k1, ns); freemem(k2, ns); freemem(k3, ns); freemem(k4, ns); freemem(k5, ns); freemem(y, ns) End; {odeiv2} End.