Program Integrals;
{$E-,N+,G+,D+}
Uses CRT;
var
  a, b, c, k : Extended;
  hx, hz, hy, y0, ye, x0, xe : Extended;
  i : Longint;
  Pi : Extended;
Function Pi_ : Extended;Assembler;
asm
  fldpi
end;
Function koeff : Extended;Assembler;
asm
  fld   a                 {a}
  fld1                    {1\a}
  fld1                    {1\1\a}
  fld1                    {1\1\1\a}
  fadd                    {2\1\a}
  fdiv  st(2), st(0)      {2\1\(a/2)}
  fadd  st(1), st(0)      {2\3\(a/2)}
  fxch                    {3\2\(a/2)}
  fdivp st(1), st(0)      {(2/3)\(a/2)}
  f2xm1                   {2^(2/3)-1\(a/2)}
  fsqrt                   {sqrt(2^(2/3)-1)\(a/2)}
  fmul  st(0), st(1)      {sqrt(2^(2/3)-1)*(a/2)}
end;
Function sqr_r(z : Extended) : Extended;
begin
  sqr_r :=  exp(   (1/3)*ln(sqr(a*z*z))   )  - z*z;
end;
Function r(z : Extended) : Extended;
begin
  r := sqrt( exp(   (1/3)*ln(sqr(a*z*z))   )  - z*z );
end;
Function z(x, y : Extended) : Extended;
begin
       z := sqr(x)+sqr(y){c*sqrt(sqr(x/a)+sqr(y/b))};
end;
Function Trapezoids(N : Longint): Extended;
var
  S, y, z : Extended;
  i, j, k : Longint;
begin
  hz := c/N;
  S := 0;
  for k := 0 to N do
  begin
   z := hz*k;
   ye := b*sqr(z/c);
   y0 := -ye;
   hy := (ye-y0)/N;
   for i := 1 to N-1 do
   begin
     y := y0+hy*i;
     xe := a*sqrt(sqr(z/c)-sqr(y/b));{}
     x0 := -xe;
     hx := (xe-x0)/N;{2*r(z0+hz*i)/N}
     for j := 1 to N-1 do
       S := S + (sqr(x0+hx*j)+sqr(y))*hy*hx*hz;{2 - две поверхности};
   end;
  end;
  Trapezoids := S;
end;
Function Simpson(N : Longint): Extended;
var
  S, y : Extended;
  i, j : Longint;
begin
  S := 0;
  hy := (ye-y0)/N;
  for i := 0 to (N div 2 - 1) do
  begin
    y := y0+hy*(2*i+1);
    xe := r(y);{}
    x0 := -xe;
    hx := (xe-x0)/N;{2*r(z0+hz*i)/N}
    for j := 0 to (N div 2 - 1) do
      S := S + 8*z(x0+(2*j+1)*hx, y)*hy*hx/3;{2 - две поверхности}
    for j := 1 to (N div 2 - 1) do
      S := S + 4*z(x0+2*j*hx, y)*hy*hx/3;{2 - две поверхности}
  end;
  for i := 1 to (N div 2 - 1) do
  begin
    y := y0+hy*2*i;
    xe := r(y);{}
    x0 := -xe;
    {Ничего не прибавляем - это нули}
    hx := (xe-x0)/N;{2*r(z0+hz*i)/N}
    for j := 0 to (N div 2 - 1) do
      S := S + 8*z(x0+(2*j+1)*hx, y)*hy*hx/3;{2 - две поверхности}
    for j := 1 to (N div 2 - 1) do
      S := S + 4*z(x0+2*j*hx, y)*hy*hx/3;{2 - две поверхности}
  end;
  Simpson := S;
end;
Function MonteCarlo(N : Longint): Extended;
var
  i : Longint;
  x, y, zz, aa, bb : Extended;
  finds : Longint;
begin
  finds := 0;
  aa := exp(1/3*ln(a));
  bb := exp(1/3*ln(b));
  for i := 1 to N do
  begin
    x := a*(2*Random-1.0);
    y := b*(2*Random-1.0);
    zz := c*Random;
    if zz>=z(x, y) then inc(finds);
  end;
  MonteCarlo := 4*a*b*c*(finds/N);
end;
Function Trapezza(N : Longint): Extended;
var
  S : Extended;
  i : Longint;
begin
  hy := (ye-y0)/N;
  S := sqr_r(ye);
    for i := 1 to N-1 do
      S := S + sqr_r(y0+i*hy)*hy;
  Trapezza := Pi*S;
end;
Function TheSimpsons(N : Longint): Extended;
var
  S : Extended;
  i : Longint;
begin
  hy := (ye-y0)/N;
  S := sqr_r(ye)*hy/3;
   for i := 0 to (N div 2-1) do
     S := S + 4*sqr_r(y0+(2*i+1)*hy)*hy/3;
   for i := 1 to (N div 2-1) do
     S := S + 2*sqr_r(y0+2*i*hy)*hy/3;
  TheSimpsons := Pi*S;
end;
begin
  ClrScr;
  Randomize;
  Pi := 0.0;
  Pi := Pi_;
  a := 2.0;
  b := 3.0;
  c := 20/78;{1/78}
  k := koeff;{a*sqrt(2^(2/3)-1)/2}
  WriteLn('k= ', k:1:20);
  WriteLn;
  WriteLn;
  WriteLn('Теоретическое значение: ', '':0, Pi*a*b*c*(sqr(a)+sqr(b))/20:1:20);
  WriteLn;
  WriteLn('Метод трапеций: ', '':8, Trapezoids(100):1:20);
  WriteLn('Метод Симпсона: ', '':8, Simpson(1000):1:20);
  WriteLn('Метод Монте-Карло: ', '':5, MonteCarlo(1000000):1:20);
  WriteLn;
  WriteLn('Метод трапецца: ', '':8, Trapezza(10000000):1:20);
  WriteLn('Метод Симпсонов: ', '':7, TheSimpsons(10000000):1:20);
  ReadKey;
end.

Назад