Метод Ньютона

[Dram]

Метод Ньютона

Попросили посмотреть программку. и сказать что она не рисует графики...посмотрел и понял что я основательно подзабыл "Паскаля" .
может подскажете что тут не так ?

uses crt,graph;
var
 gd,gm :integer;
 c1 :real;
 x,y :array [1..10] of real;
 i :longint;
 ekX,tpX :real;
 ekY,tpy :real;
 pogr :real;
 sh1,sh2,sh3 :real;
 zn1,zn2 :real;
 kolwo :longint;
 shag :real;
 int1,int2 :real;
 OsiX,osiY :integer;
 mx,my :real;
 xx,yy :integer;
 ch :string;
function ap(x:real):real;
begin
 ap:=x*c1;
end;
function l(q:real):real;
var
 i,j :integer;
 s,p :real;
begin
 s:=0;
 for i:=1 to 10 do
 begin
 p:=1;
 for j:=1 to 10 do
 begin
 if i<>j then p:=p*(q-X[j])/(X[i]-X[j]);
 end;
 s:=s+Y[i]*p
 end;
 l:=s;
end;
function razn(x:real):real;
begin
 razn:=l(x)-ap(x);
end;
function ekstr(NashPoiska,KonPoiska,e:real):real;
var
 x1,x2 :real;
 poprawka :real;
begin
 repeat
 poprawka:=(KonPoiska-NashPoiska)*0.616;
 x1:=NashPoiska+poprawka;
 x2:=KonPoiska-poprawka;
 if (l(x1)>=l(x2)) then NashPoiska:=x2 else KonPoiska:=x1;
 until abs(NashPoiska-KonPoiska)<e;
 ekstr:=NashPoiska;
end;
function kasat(a1,e:real):real;
var
 x1,x2 :real;
begin
 x2:=a1;
 repeat
 a1:=x2;
 x2:=a1-(razn(a1)*e/((razn(a1+e))-razn(a1)));
 until abs(x2-a1)<=e;
 kasat:=x2;
end;
begin
 clrscr;
 pogr:=0.01;
 writeln('Uzli interpoljacii');
 writeln(' X Y');
 for i:=1 to 10 do
 begin
 x[i]:= 2*sin((-1)*i + i*i)/4 + 2*(i-5);
 y[i]:= 2*cos((-1)*i-5) +2*sin((-1)*i+5)+ i - 5;
 writeln(x[i]:3:3,' ',y[i]:3:3);
 end;
 for i:=1 to 10 do
 begin
 sh1:=sh1+10*(X[i]*Y[i]);
 sh2:=sh2+X[i];
 sh3:=sh3+Y[i];
 zn1:=zn1+10*(X[i]*X[i]);
 zn2:=zn2+X[i];
 end;
 zn2:=zn2*zn2;
 c1:=(sh1-sh2*sh3)/(zn1-zn2);
 sh1:=0;
 sh2:=0;
 for i:=1 to 10 do
 begin
 sh1:=sh1+Y[i];
 sh2:=sh2+X[i];
 end;
 writeln('Formula approksimacii y=',c1:3:3,'*X');
 ekX:=ekstr(-5,0,pogr);
 ekY:=l(ekX);
 tpX:=kasat(-1,pogr);
 tpY:=l(tpX);
 writeln('toshka ekstremuma [',ekx:3:3,';',ekY:3:3,']');
 writeln('toshka peresheshenija [',tpx:3:3,';',tpY:3:3,']');
 kolwo:=1;
 int1:=1;
 while abs(int1-int2)>pogr*15 do
 begin
 int2:=int1;
 int1:=0;
 kolwo:=kolwo*2;
 shag:=(ekx-tpX)/kolwo;
 for i:=1 to kolwo do
 begin
 int1:=int1+shag/6*(l(tpX+shag*i)+4*l(tpX+shag*(i+0.5))+l(tpX+shag*(i+1)));
 end;
 end;
 writeln('Integral wishislenij ',int1:3:3);
 writeln('Integral wishislenij proslim progonom ',int2:3:3);
 writeln('KOlishestwo razbienij ',kolwo);
 writeln('shag razbienij ',shag:3:3);
 readkey;
 initgraph(gd,gm,'c:\bp\bgi');
 mY:=20;
 mX:=25;
 OsiX:=320;
 OsiY:=240;
 setcolor(15);
 line(0,OsiY,640,OsiY);
 line(OsiX,0,OsiX,480);
 setcolor(15);
 line(640,osiY,636,osiY-4);
 line(640,osiY,636,osiY+4);
 line(OsiX,0,osiX-3,4);
 line(OsiX,0,osiX+3,4);
 outtextXY(OsiX+5,OsiY+5,'0');
 outtextXY(OsiX+5,5,'Y');
 outtextXY(615,OsiY-15,'X');
 setcolor(15);
 for i:=1 to 10 do
 begin
 xx:=OsiX+round(X[i]*mX);
 yy:=OsiY-round(Y[i]*mY);
 circle(xx,yy,3);
 setfillstyle(1,3);
 floodfill(xx,yy,15);
 end;
 setcolor(3);
 circle(OsiX+round(tpx*mX),OsiY-round(tpy*mY),2);
 setfillstyle(1,3);
 floodfill(OsiX+round(tpx*mX),OsiY-round(tpy*mY),3);
 setcolor(5);
 circle(OsiX+round(ekx*mX),OsiY-round(eky*mY),2);
 setfillstyle(1,5);
 floodfill(OsiX+round(ekx*mX),OsiY-round(eky*mY),5);
 setlinestyle(1,1,1);
 setcolor(15);
 for i:=-20 to 20 do
 if i<>0 then
 begin
 line(osiX+round(i*mX),0,OsiX+round(i*mX),480);
 moveto(OsiX-8+round(i*2*mx+7),OsiY-15);
 setcolor(red);
 str(i*2,ch);
 outtext(ch);
 setcolor(15);
 end;
 for i:=-20 to 20 do
 if i<>0 then
 begin
 line(0,OsiY-round(i*2*mY),640,OsiY-round(i*2*mY));
 moveto(OsiX-20,OsiY-round(i*4*my));
 str(i*4,ch);
 setcolor(red);
 outtext(ch);
 setcolor(15);
 end;
 setcolor(7);
 setlinestyle(0,0,0);
 shag:=-15;
 repeat
 shag:=shag+0.001;
 yy:=round(ap(shag)*mY);
 putpixel(round(shag*mX)+OsiX,OsiY-yy,green);
 yy:=round(l(shag)*mY);
 putpixel(round(shag*mX)+OsiX,OsiY-yy,blue);

 until shag>15;
 setcolor(blue);
 line(round(tpX*mx)+OsiX,OsiY-round(tpY*my),round(tpX*mx)+OsiX,OsiY);
 line(round(ekX*mx)+OsiX,OsiY-round(ekY*my),round(ekX*mx)+OsiX,OsiY);
 line(round(tpX*mx)-1+OsiX,OsiY,round(ekX*mx)+1+OsiX,OsiY);
 setfillstyle(2,blue);
 floodfill(round(ekX*mx)+OsiX+2,OsiY+2,blue);
 readkey;
end.

 

[Vladimir_S]

Лёша, как-то не очень тянет влезать "с потрохами" в чужую программу. Хотя немного поиграл. Могу сказать, что ошибка появляется при первом же обращении к функции l(q) - скорее всего, там возникает деление на 0. Нужно вывести все значения массива Х и посмотреть, нет ли совпадающих. Если нет, двигаться дальше. Вероятнее всего, внутри всё той же функции.