mirror of
https://gitlab.com/freepascal.org/fpc/source.git
synced 2025-04-05 14:48:18 +02:00
427 lines
9.8 KiB
ObjectPascal
427 lines
9.8 KiB
ObjectPascal
program InitInfEx;
|
||
|
||
uses
|
||
Typ,
|
||
Spe,
|
||
Int;
|
||
|
||
var
|
||
num2: ArbInt;
|
||
inte: ArbFloat;
|
||
const
|
||
e = 2.71828182845905;
|
||
|
||
function cx(x: ArbFloat): ArbFloat;
|
||
begin
|
||
cx := cos(x) / (sqr(x) + 1);
|
||
end;
|
||
|
||
function pcx(x: ArbFloat): ArbFloat;
|
||
begin
|
||
pcx := 2 * x * sin(x) / sqr(sqr(x) + 1);
|
||
end;
|
||
|
||
function ppcx(x: ArbFloat): ArbFloat;
|
||
var
|
||
s: ArbFloat;
|
||
begin
|
||
s := sqr(x);
|
||
ppcx := cos(x) * (2 - 6 * s) / ((s + 1) * sqr(s + 1));
|
||
end;
|
||
|
||
function cc2(x: ArbFloat): ArbFloat;
|
||
var
|
||
x2: ArbFloat;
|
||
begin
|
||
x2 := sqr(x);
|
||
cc2 := cos(x) / (sqr(x2) + x2 + 1);
|
||
end;
|
||
|
||
function ucx(x: ArbFloat): ArbFloat;
|
||
begin
|
||
if x = 0 then
|
||
ucx := 0
|
||
else
|
||
ucx := cx((1 - x) / x) / sqr(x);
|
||
end;
|
||
|
||
function ucc2(x: ArbFloat): ArbFloat;
|
||
begin
|
||
if x = 0 then
|
||
ucc2 := 0
|
||
else
|
||
ucc2 := cc2((1 - x) / x) / sqr(x);
|
||
end;
|
||
|
||
function uz(x: ArbFloat): ArbFloat;
|
||
begin
|
||
uz := sin(x) * exp(-x);
|
||
end;
|
||
|
||
function ss1(x: ArbFloat): ArbFloat; { f(«n(n-1))=0 (n=1,2,3,...) }
|
||
var
|
||
n, s, c: ArbFloat; { f(«ný)=2/(ný(n+1)) }
|
||
begin { overigens: f linear interpoleren}
|
||
s := sqrt(2 * x);
|
||
n := trunc(s);
|
||
if n * (n + 1) / 2 <= x then
|
||
n := n + 1; { n z.d.d. «n(n-1) ó x ó «n(n-1) }
|
||
c := 4 / (n * sqr(n) * (n + 1));
|
||
if s < n then
|
||
ss1 := c * (x - n * (n - 1) / 2)
|
||
else
|
||
ss1 := c * (n * (n + 1) / 2 - x);
|
||
end;
|
||
|
||
function ss2(x: ArbFloat): ArbFloat; { als ss1 met f(«ný)=2/(n.2ü) }
|
||
var
|
||
n, s, c: ArbFloat;
|
||
begin
|
||
s := sqrt(2 * x);
|
||
n := trunc(s);
|
||
if n * (n + 1) / 2 <= x then
|
||
n := n + 1;
|
||
c := spepow(2, 2 - n) / sqr(n);
|
||
if s < n then
|
||
ss2 := c * (x - n * (n - 1) / 2)
|
||
else
|
||
ss2 := c * (n * (n + 1) / 2 - x);
|
||
end;
|
||
|
||
function ss3(x: ArbFloat): ArbFloat; { x z.d.d. «.2ü ó x+1 ó 2ü s=2ü}
|
||
var
|
||
s, c, f, x1: ArbFloat;
|
||
n: ArbInt; {n even: f(n)=-4/(n.2ü)}
|
||
begin {n oneven: f(n) = 4/(n.2ü)}
|
||
n := 0;
|
||
s := 1;
|
||
x1 := x + 1; { overigens: f lineair interpol.}
|
||
repeat
|
||
n := n + 1;
|
||
s := s * 2
|
||
until s > x1;
|
||
c := 16 / (n * sqr(s));
|
||
if x1 < 0.75 * s then
|
||
f := c * (x1 - s / 2)
|
||
else
|
||
f := c * (s - x1);
|
||
if odd(n) then
|
||
ss3 := f
|
||
else
|
||
ss3 := -f;
|
||
end;
|
||
|
||
function ss4(x: ArbFloat): ArbFloat; { 0 ó x ó 1}
|
||
var
|
||
y, h: ArbFloat; { zij x = ä [1:ì] c(n)/3ü c(n)=0,1,2}
|
||
ready: boolean; { zoek kleinste k met c(k)=1}
|
||
begin { dan geldt f(x)=u(y)/3k }
|
||
y := 3 * x;
|
||
h := 1 / 3;
|
||
ready := False; { met u(y)=|y-1«|, 1 ó y ó 2}
|
||
repeat { en y = ä [k:ì] (c(n)/3ü)*3k }
|
||
if (y < 1) or (y > 2) then
|
||
begin
|
||
if y < 1 then
|
||
y := 3 * y
|
||
else
|
||
y := 3 * (y - 2);
|
||
h := h / 3;
|
||
if h < macheps then
|
||
begin
|
||
ready := True;
|
||
ss4 := 0;
|
||
end;
|
||
end
|
||
else
|
||
begin
|
||
ready := True;
|
||
if y < 1.5 then
|
||
ss4 := h * (y - 1)
|
||
else
|
||
ss4 := h * (2 - y);
|
||
end
|
||
until ready;
|
||
end;
|
||
|
||
function ss5(x: ArbFloat): ArbFloat; { uitbreiding ss4}
|
||
var
|
||
y, h: ArbFloat; {functiewaarden op 'volgend' interval}
|
||
ready: boolean; { [n, n+1] telkens halveren}
|
||
begin
|
||
y := 3 * frac(x);
|
||
h := spepow(0.5, trunc(x)) / 3;
|
||
ready := False;
|
||
repeat
|
||
if (y < 1) or (y > 2) then
|
||
begin
|
||
if y < 1 then
|
||
y := 3 * y
|
||
else
|
||
y := 3 * (y - 2);
|
||
h := h / 3;
|
||
if h < macheps then
|
||
begin
|
||
ready := True;
|
||
ss5 := 0;
|
||
end;
|
||
end
|
||
else
|
||
begin
|
||
ready := True;
|
||
if y < 1.5 then
|
||
ss5 := h * (y - 1)
|
||
else
|
||
ss5 := h * (2 - y);
|
||
end
|
||
until ready;
|
||
end;
|
||
|
||
function ss6(x: ArbFloat): ArbFloat; { 0 ó x ó 1, 'gladdere' variant van ss4}
|
||
var
|
||
y, h: ArbFloat;
|
||
ready: boolean;
|
||
|
||
function f(y: ArbFloat): ArbFloat; { 1 ó y ó 2, 1 x cont. diff, symm. max in 1.5}
|
||
begin
|
||
if y > 1.5 then
|
||
y := 3 - y;
|
||
if y < 1.25 then
|
||
f := sqr(y - 1)
|
||
else
|
||
f := 0.125 - sqr(1.5 - y);
|
||
end;
|
||
|
||
begin
|
||
y := 3 * x;
|
||
h := 1 / 3;
|
||
ready := False;
|
||
repeat
|
||
if (y < 1) or (y > 2) then
|
||
begin
|
||
if y < 1 then
|
||
y := 3 * y
|
||
else
|
||
y := 3 * (y - 2);
|
||
h := h / 3;
|
||
if h < macheps then
|
||
begin
|
||
ready := True;
|
||
ss6 := 0;
|
||
end;
|
||
end
|
||
else
|
||
begin
|
||
ready := True;
|
||
ss6 := h * f(y);
|
||
end
|
||
until ready;
|
||
end;
|
||
|
||
function bb(x: ArbFloat): ArbFloat;
|
||
begin
|
||
bb := spepow(x, -x) * (ln(x) + 1);
|
||
end;
|
||
|
||
var
|
||
integral, ae, err: ArbFloat;
|
||
term: ArbInt;
|
||
intex: boolean;
|
||
|
||
procedure Header;
|
||
begin
|
||
Write('int': num2, '': numdig - num2, ' ', 'err': 7, ' ': 4);
|
||
if intex then
|
||
Write('f': 6);
|
||
writeln;
|
||
end;
|
||
|
||
procedure ShowResults;
|
||
var
|
||
f: ArbFloat;
|
||
begin
|
||
if intex then
|
||
f := inte - integral;
|
||
case term of
|
||
1:
|
||
begin
|
||
Write(integral: numdig, ' ', err: 10, ' ');
|
||
if intex then
|
||
writeln(f: 10)
|
||
else
|
||
writeln;
|
||
end;
|
||
2:
|
||
begin
|
||
Write(integral: numdig, ' ', err: 10, ' ');
|
||
if intex then
|
||
writeln(f: 10)
|
||
else
|
||
writeln;
|
||
Writeln(' process afgebroken, te hoge nauwkeurigheid?');
|
||
end;
|
||
3: Writeln('Verkeerde parameterwaarde (<=0) bij aanroep: ', ae: 8);
|
||
4:
|
||
begin
|
||
Write(integral: numdig, ' ', err: 10, ' ');
|
||
if intex then
|
||
writeln(f: 10)
|
||
else
|
||
writeln;
|
||
writeln(' process afgebroken, moeilijk, mogelijk divergent?');
|
||
end;
|
||
end;
|
||
end;
|
||
|
||
begin
|
||
num2 := numdig div 2;
|
||
|
||
Writeln(' ì ');
|
||
Writeln(' ô cos x ã ');
|
||
Writeln(' ³ ------- dx = -- ');
|
||
Writeln(' 0 õ xý+ 1 2e ');
|
||
writeln;
|
||
ae := 1e-8;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
writeln;
|
||
inte := 0.5 * pi / e;
|
||
intex := True;
|
||
writeln(inte: numdig, ' is ''exacte'' oplossing');
|
||
Int1fr(@cx, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
ShowResults;
|
||
writeln;
|
||
Writeln('berekend met Int1fr via transformatie x=(1-t)/t');
|
||
writeln;
|
||
Int1fr(@ucx, 0, 1, ae, integral, err, term);
|
||
Header;
|
||
ShowResults;
|
||
Writeln(' ì ');
|
||
Writeln(' ô 2x sin x ã ');
|
||
Writeln(' ³ --------- dx = -- ');
|
||
Writeln(' 0 õ (xý+ 1)ý 2e ');
|
||
writeln;
|
||
ae := 1e-8;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
writeln;
|
||
Int1fr(@pcx, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
ShowResults;
|
||
Writeln(' ì ');
|
||
Writeln(' ô (2-6xý)cos x ã ');
|
||
Writeln(' ³ ------------ dx = -- ');
|
||
Writeln(' 0 õ (xý+ 1)3 2e ');
|
||
writeln;
|
||
ae := 1e-8;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
writeln;
|
||
Int1fr(@ppcx, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
ShowResults;
|
||
|
||
Writeln(' ì ');
|
||
Writeln(' ô cos x ');
|
||
Writeln(' ³ ------------ dx = (ã/û3) exp(-«û3) sin(ã/6+«) ');
|
||
Writeln(' 0 õ (xý)ý+xý+ 1 ');
|
||
writeln;
|
||
writeln;
|
||
ae := 1e-8;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
writeln;
|
||
inte := (pi / sqrt(3)) * exp(-sqrt(0.75)) * sin(pi / 6 + 0.5);
|
||
writeln(inte: numdig, ' is ''exacte'' oplossing');
|
||
Int1fr(@cc2, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
ShowResults;
|
||
writeln;
|
||
Writeln('berekend met Int1fr via transformatie x=(1-t)/t');
|
||
writeln;
|
||
writeln;
|
||
writeln(inte: numdig, ' is ''exacte'' oplossing');
|
||
Int1fr(@ucc2, 0, 1, ae, integral, err, term);
|
||
Header;
|
||
ShowResults;
|
||
|
||
Writeln(' ì ');
|
||
Writeln(' ô sin u ');
|
||
Writeln(' ³ ------ du = « ');
|
||
Writeln(' õ exp(u) ');
|
||
writeln(' 0 ');
|
||
ae := 1e-8;
|
||
intex := True;
|
||
inte := 0.5;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
writeln;
|
||
Int1fr(@uz, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
ShowResults;
|
||
|
||
writeln(' functie ss1; int = ä {1:ì}1/n(n+1) = 1');
|
||
ae := 1e-8;
|
||
intex := True;
|
||
inte := 1;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
int1fr(@ss1, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
Showresults;
|
||
|
||
writeln(' functie ss2; int = ä {1:ì} («)ü = 1');
|
||
ae := 1e-8;
|
||
intex := True;
|
||
inte := 1;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
int1fr(@ss2, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
Showresults;
|
||
|
||
writeln(' functie ss3; int = ä {1:ì} (-1)ü/n = ln(2)');
|
||
ae := 1e-8;
|
||
intex := True;
|
||
inte := ln(2);
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
int1fr(@ss3, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
Showresults;
|
||
|
||
writeln(' functie ss4 (op [0,1]); int = 1/28 ');
|
||
ae := 1e-8;
|
||
intex := True;
|
||
inte := 1 / 28;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
int1fr(@ss4, 0, 1, ae, integral, err, term);
|
||
Header;
|
||
Showresults;
|
||
|
||
writeln(' functie ss5 (op [0,ì)); int = 1/14 ');
|
||
ae := 1e-8;
|
||
intex := True;
|
||
inte := 1 / 14;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
int1fr(@ss5, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
Showresults;
|
||
|
||
writeln(' functie ss6 (op [0,1]); int = 1/112 ');
|
||
ae := 1e-8;
|
||
intex := True;
|
||
inte := 1 / 112;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
int1fr(@ss6, 0, 1, ae, integral, err, term);
|
||
Header;
|
||
Showresults;
|
||
|
||
Writeln(' ì ');
|
||
Writeln(' ô ln(x)+1 ');
|
||
Writeln(' ³ ---------- dx = 1 ');
|
||
Writeln(' õ xx ');
|
||
writeln(' 1 ');
|
||
ae := 1e-8;
|
||
intex := True;
|
||
inte := 1;
|
||
Writeln(' Gevraagde nauwkeurigheid', ae: 12);
|
||
writeln;
|
||
Int1fr(@bb, 0, infinity, ae, integral, err, term);
|
||
Header;
|
||
ShowResults;
|
||
end.
|