fpc/packages/numlib/tests/odeiv2te.pas
2016-03-05 19:16:54 +00:00

77 lines
1.6 KiB
ObjectPascal

program odeiv2te;
uses
typ,
ode;
const
n1 = 3;
n2 = 4;
n3 = 6;
n = n2 - n1 + 1;
n4 = n3 + n - 1;
var
ex, nv, i, j, k, h, term: ArbInt;
a, b, d, ae: ArbFloat;
ya: array[n1..n2] of ArbFloat;
yb: array[n3..n4] of ArbFloat;
procedure f(x: ArbFloat; var y, y1: ArbFloat);
var
yloc: array[1..n] of ArbFloat absolute y;
y1loc: array[1..n] of ArbFloat absolute y1;
begin
y1loc[1] := 2 * x * yloc[1] + yloc[2];
y1loc[2] := -yloc[1] + 2 * x * yloc[2];
end; {f}
function phi1(x: ArbFloat): ArbFloat;
begin
phi1 := exp(x * x) * sin(x);
end; {phi1}
function phi2(x: ArbFloat): ArbFloat;
begin
phi2 := exp(x * x) * cos(x);
end; {phi2}
begin
Write(' program results odeiv2te');
case sizeof(ArbFloat) of
4: writeln('(single)');
6: writeln('(real)');
8: writeln('(double)');
end;
Read(nv);
writeln;
writeln(' number of examples: ', nv: 2);
k := numdig;
h := k div 2;
for ex := 1 to nv do
begin
writeln;
writeln(' example number :', ex: 2);
Read(a, b);
for j := n1 to n2 do
Read(ya[j]);
Read(ae);
d := b - a;
writeln;
writeln(' ae =', ae: 10);
writeln;
writeln('b': 3, 'yb1': h + 4, 'yb2': k, 'phi1(x)': k + 3, 'phi2(x)': k - 2, 'term': h + 2);
for i := 1 to 5 do
begin
odeiv2(@f, a, ya[n1], b, yb[n3], n, ae, term);
writeln(b: 5: 2, yb[n3]: k, yb[n3 + 1]: k, phi1(b): k, phi2(b): k, term: 3);
a := b;
for j := n1 to n2 do
ya[j] := yb[n3 - n1 + j];
b := b + d;
end; {i}
writeln(' -------------------------------------------------');
end; {ex}
Close(input);
Close(output);
end.