mirror of
https://gitlab.com/freepascal.org/fpc/source.git
synced 2025-07-15 23:46:03 +02:00
342 lines
9.4 KiB
ObjectPascal
342 lines
9.4 KiB
ObjectPascal
{
|
|
$Id$
|
|
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.
|
|
{
|
|
$Log$
|
|
Revision 1.2 2000-01-25 20:21:42 marco
|
|
* small updates, crlf fix, and RTE 207 problem
|
|
|
|
Revision 1.1 2000/01/24 22:08:58 marco
|
|
* initial version
|
|
|
|
|
|
}
|