Код:
uses
wincrt;
const
mu = 0.; sig = 1.;
eps = 0.1e-3; { точность решения уравнения F((ln(x) - mk) / (sig * sqrt(2))) = p }
smax = 0.7; sk = 0.2;
kh0 = 0.7; kb0 = 0.5;
type
func= function(x: extended): extended;
var
ft: text;
name: string;
mk, dk, s, h, p, kp: extended;
jh, jb, khf, kbf: extended;
{ функции Ошибок }
function erf(x: extended): extended;
const
eps = 0.1e-3; { точность вычисления функции Ошибок }
var
sn1, sn2: extended;
n: integer;
{ функция rowerf, вызываемая в функции erf, }
{ для вычисления n членов ряда функции Ошибок }
function rowerf(n: integer; x: extended): extended; far;
var
a, b, c, s: extended;
i: integer;
begin
a:=x; b:=1; c:=1; s:=x;
for i:=1 to n do
begin
a:=-a * sqr(x); b:=b + 2; c:=c * i;
s:=s + a / (b * c);
end;
rowerf:= 2 / sqrt(pi) * s
end; { of rowerf }
begin
n:=8;
repeat
sn1:=rowerf(n, x); sn2:=rowerf(n + 1, x);
inc(n)
until (abs(sn2 - sn1) <= eps);
erf:=sn2
end; { of erf }
{ функция вычисления определенного интеграла по методу Симпсона }
function integs(f: func; a,b: extended): extended; far;
const
eps = 0.001; { точность вычисления определенного интеграла }
var
s, h, x: extended;
ds1, ds2, ds, int: extended;
n: integer;
begin
n:=10; ds1:=0;
repeat
s:=0; h:=(b - a) / n; x:=a;
repeat
s:=s + (f(x) + 4 * f(x + h * 0.5) + f(x + h)) / 6;
x:=x + h
until (x > b - h);
int:=h * s;
ds2:=int; ds:=abs(ds2 - ds1); ds1:=ds2;
n:=2 * n
until (ds <= eps);
integs:=int
end; { of integs }
{ функция распределения логарифмически-нормального закона }
function Ferf(x: extended): extended; far;
begin
Ferf:=(1 + erf(x)) * 0.5
end; { of Ferf }
{ функция решения уравнения F((ln(x) - mk) / (sig * sqrt(2))) = p }
function fkp(Ferf: func; p: extended): extended;
var
x, r: extended;
begin
x:=eps; r:=(ln(x) - mk) / (sig * sqrt(2));
while (abs(Ferf(r) - p) > eps) do
begin
x:=x + eps;
r:=(ln(x) - mk) / (sig * sqrt(2))
end;
fkp:=x
end; { of fkp }
{ x * f(x) (f(x) - функция плотности распределения логарифмически-нормального закона) }
function f(x: extended): extended; far; { подинтегральная функция }
begin
f:=x / (sqrt(2 * pi) * sig * x) * (exp(-sqr(ln(x) - mk) / 2 * sqr(sig)))
end; { of f }
begin
write(' Введите полное имя файла результата: '); readln(name);
assign(ft, name);
rewrite(ft);
writeln(' Относительные фазовые проницаемости лог-нормального закона распределения для нефтенасыщенности');
writeln;
writeln(ft,' Относительные фазовые проницаемости лог-нормального закона распределения для нефтенасыщенности');
writeln(ft);
{writeln(' mu = ', mu:6:3,' sig = ', sig:6:3); readln;
writeln(ft, ' mu = ', mu:6:3,' sig = ', sig:6:3);
mk:=exp(mu + sqr(sig) / 2);}
mk:=0.255;
mk:=0.3;
dk:=sqr(mk) * (exp(sqr(sig)) - 1);
{ mk - математическое ожидание, dk - дисперсия }
writeln(' mk = ', mk:6:3,' dk = ', dk:6:3);
writeln; readln;
writeln(ft, ' mk = ', mk:6:3,' dk = ', dk:6:3);
writeln(ft);
s:=sk; h:=0.1;
repeat
writeln(' s = ', s:6:3);
writeln(ft, ' s = ', s:6:3);
p:=(s - sk) / (smax - sk);
if (p = 0) then begin
jh:=0;
jb:=mk
end
else begin
kp:=fkp(Ferf, p);
writeln(' kp = ', kp:6:3);
writeln(ft, ' kp = ', kp:6:3);
writeln(' p =', p:6:3,' Ferf = ', Ferf((ln(kp) - mk) / (sig * sqrt(2))):6:3);
writeln(ft, ' p =', p:6:3,' Ferf = ', Ferf((ln(kp) - mk) / (sig * sqrt(2))):6:3);
jh:=integs(f, eps, kp);
if (jh > mk) then jh:=mk;
jb:=mk - jh
end;
writeln(' jh = ',jh:6:3,' jb = ',jb:6:3);
writeln(ft, ' jh = ',jh:6:3,' jb = ',jb:6:3);
khf:=kh0 * (jh / mk); kbf:=kb0 * (jb / mk);
writeln(' khf = ', khf:6:3,' kbf = ', kbf:6:3);
writeln(ft, ' khf = ', khf:6:3,' kbf = ', kbf:6:3);
writeln; writeln(ft); readln;
s:=s + h
until (s > smax - h / 2);
writeln(' s = ',s:6:3);
writeln(ft, ' s = ',s:6:3);
jh:=mk;
jb:=0;
writeln(' jh = ',jh:6:3,' jb = ',jb:6:3);
writeln(ft, ' jh = ',jh:6:3,' jb = ',jb:6:3);
khf:=kh0 * (jh / mk); kbf:=kb0 * (jb / mk);
writeln(' khf = ',khf:6:3,' kbf = ',kbf:6:3);
writeln(ft, ' khf = ',khf:6:3,' kbf = ',kbf:6:3);
close(ft);
writeln(' Результат см. в файле ', name)
end.