Численные методы

uses graph;
var gd,gm,i,n:integer; x,y,c1,c2,xn,yn,x2n,xyn:real;  k,m: array [1..100] of real;
  function fx(x:real):integer;
begin
     fx:=round(320+150*x);
end;
  function fy(y:real):integer;
begin
     fy:=round(240-150*y);
end;
   function f(x:real):real;
            function P(i:integer;x:real):real;
            var j:integer; res:real;
         begin
            res:=1;
            for j:=1 to n do
            if (i<>j) then res:=res*(x-k[j])/(k[i]-k[j]);
            P:=res;
         end;
             function L(x:real):real;
             var i:integer; res:real;
         begin
              res:=0;
              for i:=1 to n do
              res:= res+m[i]*P(i,x);
              L:=res;
         end;
 
   begin
     f:=L(x);
   end;
 
procedure osi;
     begin
          line(fx(0),fy(-40),fx (0), fy(40));
          line(fx(-32),fy(0),fx (32), fy(0));
        for i:=-100 to 100 do
            line(fx(i),fy(0.2),fx (i), fy(-0.2));
        for i:=-40 to 40 do
            line(fx(-0.3),fy(i),fx (0.3), fy(i));
     end;
 
begin
     gd:=detect;
     initgraph(gd,gm,' ');
     if graphresult <> grok then halt(1);
 
n:=10;
for i:=1 to n do
    begin
         k[i]:=3*sin((-2)*i + i*i)/4 + (-2)*(i-5);
         m[i]:=3*cos((-2)*i-5) +(-2)*sin((-2)*i+5)+ i - 5;
         setcolor(red);
         circle(fx(k[i]),fy(m[i]),3);
    end;
 
for i:=1 to n do
       begin
            xn:=xn+k[i];
            xyn:=xyn+k[i]*m[i];
            yn:=yn+m[i];
            x2n:=x2n+k[i]*k[i];
       end;
xn:=xn/5;
yn:=yn/5;
xyn:=xyn/5;
x2n:=x2n/5;
c1:=(xn*yn+xyn)/(x2n-xn*xn);
c2:=yn-c1*xn;
setcolor(blue);
line(fx(-1),fy(c1*(-1)+c2),fx(1),fy(c1+c2));
 
 
     x:=-2; setcolor(white); osi;
     while (x<2) do
           begin
                putpixel(fx(x),fy(f(x)),77);
                x:=x+0.001;
 
           end;
 
      readln;
      closegraph;
end.

program prostie iteracii;
var i:integer; a,b,e:real; x:Array[0..1000] of real;
 
function u(x:real):real;
begin
u:=2-sin(1/x); {formula vida 'x = ...'}
end;
 
begin
a:=1.2; b:=2; e:=0.001; {na4alo, konec, to4nost}
i:=0;
x:=(a+b)/2; {na4alnoe priblijenie}
repeat
inc(i);
x:=u(x[i-1]); {formula}
until (abs(x-x[i-1]) <= e); {po opredeleniy}
writeln(x[i-1]:0:10); {rezultat}
 
readln; readln;
end.
 
program metodnytona;
var a,b,e:real;
 
function Fx(x:real):real; {s4itaet zna4enie fuini v to4ke x}
begin
Fx:=0.1*x*x-x*ln(x); {funkci9 f(x)}
end;
 
function Fx1(x:real):real; {s4itaet priblijennoe zna4enie proizvodnoi funkcii}
const z = 0.0000001;
begin
Fx1:=(Fx(x+z)-Fx(x))/z; {iz matana}
end;
 
procedure Newton(x0:real);
var i:integer; x:array[0..1000] of real;
begin
i:=-1;
x[0]:=x0;
repeat
inc(i);
x[i+1]:=x-Fx(x)/Fx1(x); {formula}
until (abs(x[i+1]-x) <= e); {po opredeleniy}
writeln('result: x',i,'= ',x:0:10); {rezultat}
end;
 
begin
a:=1; b:=2; e:=0.00001; {na4alo, konec, to4nost}
Newton((a+b)/2); {vizivaem proceduru}
 
readln; readln;
end.

program ekstremum;
const z= 1.618;
var a,b,e,x,x1,x2,y1,y2:real;
 
function fx(x:real):real; {s4itaet zna4enie fuini v to4ke x}
begin
fx:=sin(x/4)*cos(2*x)*(-5*x*x+2*x+6); {funkci9 f(x)}
end;
 
begin
write('a = '); readln(a); {na4alo}
write('b = '); readln(b); {konec}
write('e = '); readln(e); {to4nost}
x1:=b-(b-a)/z; x2:=a+(b-a)/z; {dve simmetri4nie to4ki otnositelno centra otrezka}
y1:=fx(x1); y2:=fx(x2); {zna4enie funkcii v etih to4kah}
repeat
{otbrasivaem tot iz koncov otrezka, k kotoromu sredi dvuh to4ek okazalas ta, zna4enie kotoroi maksimalno}
if(y1<=y2) then
begin
b:=x2;
x2:=x1;
x1:=b-(b-a)/z;
y2:=y1;
y1:=fx(x1);
end else
begin
a:=x1;
x1:=x2;
x2:=a+(b-a)/z;
y1:=y2;
y2:=fx(x2);
end;
until (abs(b-a)<=e); {povtor9em poka ne dostignem zadanoi to5nosti}
x:=(b+a)/2; {to4ka ekstremuma}
writeln('f(x)=',fx(x):0:10,' x=',x:0:10); {rezultat}
 
readln; readln;
end.

uses crt,graph;
 
{nastrojki grafika}
    width = 640;
    height = 640;
    shag_setki_x = 20;
    shag_setki_y = 1;
    x_offset = 62;
    y_offset = 407;
    masshtab_x = 2.5977653631285;
    masshtab_y = 33.990382025225;
{tocnost}
    eps = 0.01;
 
var
{uzlovie znacenija}
    x,y: array[1..10] of double;
{s4et4ik dlya ispolzovanija v raznih ciklah}
    i: integer;
{dlya graphiciskih funkcij}
    gd,gm: integer;
{krajnie zna4enija}
    min_x,max_x: double;
{dlya risovanija setki. tekushie znacenija}
    px,py,valu: integer;
    s: string;
{parametri approximacii}
    C,a_xy,a_x2: double;
{shag risovanija interpolyacii}
    sx,sy: double;
{proizv - proizvodnaya. cross_x i cross_x0 - priblizheniya
 k pereseceniyu interpolyacii i approximacii.
 hord_c - dlya perevorota promezhutka}
    proizv, cross_x, cross_x0, hord_c: double;
{dlya poiska extremuma}
    ext_a, ext_b, ext_x, ext_x1, ext_y: double;
{dlya integrirovaniya}
    integr_a, integr_b, integr: double;
    integr_n: integer;
10
{funkcii preobrazovanija obi4nih koordinat v koordinati na ekrane}
function cx(x:double):integer;
begin
    cx:=x_offset+round(masshtab_x*x);
end;
function cy(y:double):integer;
begin
    cy:=y_offset-round(masshtab_y*y);
end;
 
{approximacia}
function approx(qx:double):double;
begin
    approx:=C*qx;
end;
 
{interpoljacionnij mnogo4len lagranga}
function lagrange(qx:double):double;
var
    j:integer;
    La,l:double;
begin
    La:=0;
    for j:=1 to 10 do begin
        l:=1;
        for i:=1 to 10 do begin
            if i<>j then
                l:=l*(qx-x[i])/(x[j]-x[i]);
        end;
        La:=La+y[j]*l;
    end;
    lagrange:=La;
end;
 
{funkciya raznosti mezhdu grafikom interpolyacii i approximacii}
function diff(x:double):double;
begin
    diff:=lagrange(x)-approx(x);
end;
 
{funkcija integrirovaniya}
function integral(n:integer):double;
 
 
11
var
    integ, h: double; {zna4enie integrala i shag integrirovaniya}
    i: integer;
 
begin
    integ:=0;
    h:=(integr_b-integr_a)/n;
    {
     Xi = integr_a+i*h
     integr_a - pervaja granica intervala, i*h - itiy shag
    }
    {simpson}
    integ:=lagrange(integr_a)+lagrange(integr_b);
    for i:=1 to n do
        integ:=4*lagrange(integr_a+(i-1/2)*h)+integ;
    for i:=1 to n-1 do
        integ:=2*lagrange(integr_a+i*h)+integ;
    integ:=integ/6;
    integral:=integ*h; {umnozhaem na shag}
end;
 
begin
    clrscr;
    writeln('Uzlovie znacenija:');
    for i:=1 to 10 do begin
        x[i]:= 3*sin((-2)*i + i*i)/4 + (-2)*(i-5);
        y[i]:= 3*cos((-2)*i-5) +(-2)*sin((-2)*i+5)+ i - 5;
        writeln('  ',x[i],' ',y[i]);
    end;
    readln;
 
    {
     scitaem parametri approximacii
     peremennie a_* -- eto algebraiceskie summi
    }
    a_xy:=0; a_x2:=0;
    for i:=1 to 10 do begin
       
 
12
        a_xy := a_xy + x[i]*y[i];
        a_x2 := a_x2 + x[i]*x[i];
    end;
    C:=a_xy/a_x2;
    B:=
    {
     ishem tocku pereseciniya grafika
     interpolyacii i pryamoj approximacii
    }
    {promezhutok poiska}
    cross_x0:=60;
    cross_x :=70;
    {hordi}
    {dvojnaya proizvodnaya}
    proizv:=diff(cross_x0+2*eps)-2*diff(cross_x0+eps)+diff(cross_x0);
    proizv:=proizv/(2*eps*eps);
    if proizv*diff(cross_x0)>0 then
        hord_c:=cross_x0
    else begin
        hord_c:=cross_x;
        cross_x:=cross_x0;
    end;
    while abs(cross_x-cross_x0)>eps do begin
        cross_x0:=cross_x;
        cross_x:=cross_x-diff(cross_x)*(hord_c-cross_x)/(diff(hord_c)-diff(cross_x));
    end;
    { nahodim extremum}
        ext_a:=110;
    ext_b:=120;
 
    ext_x:=ext_a;
    ext_y:=lagrange(ext_x);
    while ext_a<ext_b do begin
        ext_a:=ext_a+eps;
        sy:=lagrange(ext_a);
        if sy<ext_y then begin
            ext_x:=ext_a;
            ext_y:=sy;
        end;
    end;
{
 integrirovanie funkcii lagranga mezhdu
 tochkoj extremuma i tockoj pereseceniya
 s approximaciej
}
13
 if ext_x<cross_x then begin
        integr_a:=ext_x;
        integr_b:=cross_x;
    end else begin
        integr_a:=cross_x;
        integr_b:=ext_x;
    end;
 
    {
     OTRISOVKA
    }
    clrscr;
    gd:=Detect;
    gm:=2; {640x480}
    {inizializiruem graficeskij adapter}
    InitGraph(gd,gm,'C:\TP\BGI\');
    SetBkColor(white);
 
    {
     krajnie znacenija grafika
    }
    max_x:=shag_setki_x; min_x:=-shag_setki_x;
    for i:=1 to 10 do begin
        if max_x<x[i] then max_x:=x[i]
        else if min_x>x[i] then min_x:=x[i];
    end;
 
    {zapolnyaem promezhutok integrirovaniya}
    SetFillStyle(LtSlashFill, LightGray);
    sx:=integr_a;
    while sx<integr_b do begin
        Bar(cx(sx), cy(lagrange(sx)), cx(sx)+1, cy(0));
        sx:=sx+1/masshtab_x;
    end;
 
    {
     ustanavlivaem cvet
     dlya osej coordinat
    }
    SetColor(blue);
    {
     risuem osi coordinat
    }
    {os OX}
14
  Line(2,y_offset,width-2,y_offset);
      {strelocka}
      Line(width-2,y_offset,width-10,y_offset+2);
      Line(width-2,y_offset,width-10,y_offset-2);
    {os OY verhnyaya cast}
    Line(x_offset,2,x_offset,height-2);
      {strelocka}
      Line(x_offset,2,x_offset-2,10);
      Line(x_offset,2,x_offset+2,10);
 
{nasecki}
    SetColor(blue);
    {po X ot 0 do +inf}
    px:=cx(shag_setki_x);
    valu:=shag_setki_x;
    while px<(width-5) do begin
        Line(px,y_offset+2,px,y_offset-2);
        Str(valu,s);
        OutTextXY(px-5, y_offset+5, s);
 
        valu:=valu+shag_setki_x;
        px:=cx(valu);
    end;
    {po Y ot 0 do +inf}
    py:=cy(shag_setki_y);
    valu:=shag_setki_y;
    while py>5 do begin
        Line(x_offset+2,py,x_offset-2,py);
        Str(valu,s);
        OutTextXY(x_offset+5, py-4, s);
 
        valu:=valu+shag_setki_y;
        py:=cy(valu);
    end;
    {po X ot 0 do -inf}
    px:=cx(-shag_setki_x);
    valu:=-shag_setki_x;
    while px>5 do begin
        Line(px,y_offset+2,px,y_offset-2);
        Str(valu,s);
        OutTextXY(px-5,y_offset+5,s);
 
       
 
15
        valu:=valu-shag_setki_x;
        px:=cx(valu);
    end;
    {po Y ot 0 do -inf}
    py:=cy(-shag_setki_y);
    valu:=-shag_setki_y;
    while py<(height-5) do begin
        Line(x_offset+2,py,x_offset-2,py);
        Str(valu,s);
        OutTextXY(x_offset+5, py-4, s);
 
        valu:=valu-shag_setki_y;
        py:=cy(valu);
    end;
 
    {
     vivodim vse tocki
    }
    SetColor(magenta);
    setfillstyle(SolidFill, magenta);
    for i:=1 to 10 do begin
        FillEllipse(cx(x[i]), cy(y[i]), 2, 2);
    end;
 
    {
     ustanavlivaem sisnij cvet i
     risuem im appromixaciyu
    }
    SetColor(blue);
    Line(cx(min_x), cy(approx(min_x)), cx(max_x), cy(approx(max_x)));
 
    {
     risuem interpolaciyu s zadannoj tocnostiyu EPS
     ot pervoj do poslednej tocki.
     Kazhduyu tocku (x0+EPS, lagrange(x0+EPS)) risuem otdelno
     po tockam zelenim cvetom
    }
    sx:=min_x;
    while sx<max_x do begin
        sy:=lagrange(sx);
        PutPixel(cx(sx), cy(sy), blue);
        sx:=sx+eps;
    end;
 
16
         SetColor(blue);
    SetFillStyle(SolidFill, blue);
    FillEllipse(cx(cross_x), cy(lagrange(cross_x)), 3,3);
    FillEllipse(cx(ext_x), cy(lagrange(ext_x)), 3,3);
 
    SetColor(DarkGray);
    Line(cx(cross_x), cy(lagrange(cross_x)), cx(cross_x), cy(0));
    Line(cx(ext_x), cy(ext_y), cx(ext_x), cy(0));
 
    {povishaem cislo razbienij, poka ne dostignem otnositelnoj tocnosti EPS}
    integr_n:=4;
    repeat
        integr_n:=integr_n*2; {povishaem chislo razbienij, chtoby dobit'sa nuzhnoj tochnosti}
        integr:=integral(integr_n); {znachenie integrala}
    until abs(integr-integral(integr_n*2))<eps;
    OutTextXY(200, 10, 'Integralnaya summa interpolyacii');
    Str(integr,s);
    OutTextXY(200, 25, 'ravna: '+s);
 
    readln;
end.