Помогите с методом ньютона Опции

[Гил Бейтс]

Дали задание построить сопло и найти число маха газа в этом сопле двумя способами:методом половинного деления и мет. Ньютона.
Я записал исходные данные координат точек образования сопла, построил его, и нашел число маха M методом половинного деления, а вот Ньютоном не получается. Я так понял надо вместо вычисления точек гравика, найти производную.

program nozzle;
uses mygraph,solver1;
var
    i,j,N                           : word;
    x,y,M                           : ArrayType;
    yB                              : ArrayBigType;
    x0,x1,x2,y0,y1,y2,Gamma,Tmp,f,f1: real;
    Key                             : byte;
    c                               :real;
begin
  N:=30;
  Gamma:=1.4;
  x0:=0;x1:=1;x2:=3;
  y0:=1.5;y1:=1;y2:=2;
  for i:=1 to N do
     begin
         X[i]    := X0 + (X2-X0)*(i-1)/(n-1);
         if x[i]<x1
            then
                y[i]  := y0+(y1-y0)/(x1-x0)*(x[i]-x0)
            else
                y[i]  := y1+(y2-y1)/(x2-x1)*(x[i]-x1);
     end;
  Graphic (n,
            40,1,600,450,
            X,Y,
            11,11,
            2,2,
            'x','y','y(x)');
            f1:=y1*y1*pi;
  for i:=1 to N do
     begin
        f:=pi*sqr(y[i]);
        Tmp:=f1/f;
        if x[i]<x1
           then key:=0
           else key:=1;
         SolvEqv( Key,N,
                  Tmp,
                  Gamma,
                  M[i]);
     end;
   Graphic (n,
            40,1,600,450,
            X,M,
            11,11,
            2,2,
            'x','M','M(x)');
  for i:=1 to N do
     begin
        c := 1+(gamma-1)/2*M[i]*M[i];
        Yb[1,i] := exp(-gamma/(gamma-1)*ln©);
        Yb[2,i] := exp(-1.0/(gamma-1)*ln©);
        Yb[3,i] := exp(-1.0/2.0*ln©);
        Yb[4,i] := exp(-1.0*ln©);
     end;
   Graphic1 (4,n,
            40,1,600,450,
            M,YB,
            11,11,
            'x','M','M(x)');

end.


unit Solver1;
interface
uses crt,mygraph;
   procedure SolvEqv( Key,N   : byte;
                      Tmp     : real;
                      Gamma   : real;
                  var M       : real);
implementation
{
------------
}
Function Fun(Tmp,Gamma,M:real ):real;
var
   a,b,c : real;
begin
   a:=-(gamma+1)/(2*(gamma-1));
   b:=2/(gamma+1);
   c:=1+(gamma-1)/2*M*M;
   Fun := Tmp - exp(a*ln(b))*M*exp(a*ln©);




end;

procedure SolvEqv;
 var
    ml,mr,fl,fr,fc,mc,e:real;
    xt,ft              : ArrayType;
    i                  : word;
begin
    e:=1e-5;
   if key=0 then
    begin
      ml:=0;
      mr:=0.999;
    end
   else
    begin
      ml:=1.001;
      mr:=100;
    end;
  Fl := Fun(Tmp,Gamma,Ml);
  FR := Fun(Tmp,Gamma,MR);

  if FL*FR < 0
     then
        begin
           {writeln('ў ¤ ­­®¬ ¤Ё Ї §®­Ґ Є®аҐ­м Ґбвм');}
           {
           for i:=1 to N do
              begin
                 Xt[i] := ML + (MR-ML)*(i-1)/(n-1);
                 Ft[i] := Fun(Tmp,Gamma,Xt[i]);
              end;
           Graphic (n,
            40,1,600,450,
            Xt,ft,
            11,11,
            2,2,
            'M','F','F(M)');
            }

           Repeat
              mc:=(ml+mr)/2;
              fc:= Fun(Tmp,Gamma,MC);
              if fl*fC<0
                 then
                    begin
                       mr:=mc;
                       fr:=fC;
                       {ml:=ml;
                       fl:=fL;}
                    end
                 else
                    begin
                       ml:=mc;
                       fl:=fc;
                       {mr:=mr;
                       fr:=fr;}
                    end;
           Until (mr-ml <= e );
           M:=(mr+ml)/2;
        end
     else
           writeln('ў ¤ ­­®¬ ¤Ё Ї §®­Ґ Є®а­п ­Ґв');
end;


unit MyGraph;
Interface
uses Graph,crt;
const
    NMax = 300;
    NGrMax = 1;
type
    arraytype2=array[1..NMax] of real;
    ArrayType=array[1..NMax] of real;
    ArrayBigType=array[1..NGrMax,1..NMax] of real;
procedure Graphic(n                  : word;
                  XP1,YP1,XP2,YP2    : word;
                  X,Y                : arraytype;
                  ngr1,ngr2          : word;
                  n0x,n0y            : word;
                  namex,namey,namegr : string );
procedure Graphic1 (NGr,n              : word;
                    XP1,YP1,XP2,YP2    : word;
                    X                  : arraytype;
                    Y                  : arrayBigtype;
                    ngr1,ngr2          : word;
                    namex,namey,namegr : string );

Implementation
Procedure Graphic;
var
   XStr,YStr                        : string;
   X1,Y1,X2,Y2,Xtmp,Ytmp            : word;
   xs,ys,xf,yf                      :word;
   XMin,XMax,YMin,YMax,Xt,Yt        : real;
   GraphDriver, GraphMode, i        : Integer;
begin
   GraphDriver := Detect;
   InitGraph(GraphDriver, GraphMode, ' ');
   if GraphResult<> grOk then
     Halt(1);
   x1 := xp1 + 70;
   x2 := xp2 - 50;
   y1 := yp1 + 30;
   y2 := yp2 - 30;
   Rectangle(XP1, YP1, XP2, YP2);
   Setlinestyle(1,0,0);
   Rectangle(X1, Y1, X2, Y2);
   Outtextxy(x2+10,y2-5,namex);
   Outtextxy(x1+10,y1-15,namey);
   Outtextxy(round((x1+x2)/2),y1-15,namegr);
   XMin:=X[1];
   XMax:=X[N];
   YMin:=10000;
   YMax:=-10000;
   for i:= 1 to N do
      begin
         if(YMin>Y[i]) then YMin:=Y[i];
         if(YMax<Y[i]) then YMax:=Y[i];
          if(XMin>X[i]) then XMin:=X[i];
         if(XMax<X[i]) then XMax:=X[i];
      end;
   for i:=1 to ngr1 do
      begin
         Xtmp:= X1 + round((X2-X1)*(i-1)/(ngr1-1));
         Xt  := XMin + (XMax-XMin)*(i-1)/(ngr1-1);
         Str(Xt:6:n0x,XStr);
         OutTextXY(XTmp-35,y2+10,XStr);
         line(Xtmp,y1,Xtmp,y2);
      end;
   for i:=1 to ngr2 do
      begin
         Ytmp:=y1 + round((y2-y1)*(i-1)/(ngr2-1));
         Yt:=YMax - (YMax-YMin)*(i-1)/(ngr2-1);
         Str(Yt:6:n0y,YStr);
         OutTextXY(x1-55,YTmp,YStr);
         line(X1,Ytmp,X2,Ytmp);
      end;
   Setlinestyle(0,1,1);
   SetColor(lightcyan);
   xs := x1 + trunc((x[1]-xmin)/(xmax-xmin)*(x2-x1));
   ys := y2 - trunc((y[1]-ymin)/(ymax-ymin)*(y2-y1));
   for i:=2 to n do
      begin
         xf := x1 + trunc((x[i]-xmin)/(xmax-xmin)*(x2-x1));
         yf := y2 - trunc((y[i]-ymin)/(ymax-ymin)*(y2-y1));
         line(xs,ys,xf,yf);
         xs:= xf;
         ys := yf;
      end;
   readln;
   CloseGraph;
end;
{
------------------------
}
Procedure Graphic1;
var
   XStr,YStr                        : string;
   X1,Y1,X2,Y2,Xtmp,Ytmp            : word;
   xs,ys,xf,yf                      : word;
   XMin,XMax,YMin,YMax,Xt,Yt        : real;
   GraphDriver, GraphMode, i,j      : Integer;
begin
   GraphDriver := Detect;
   InitGraph(GraphDriver, GraphMode, ' ');
   if GraphResult<> grOk then
     Halt(1);
   x1 := xp1 + 70;
   x2 := xp2 - 50;
   y1 := yp1 + 30;
   y2 := yp2 - 30;
   Rectangle(XP1, YP1, XP2, YP2);
   Setlinestyle(1,0,0);
   Rectangle(X1, Y1, X2, Y2);
   Outtextxy(x2+10,y2-5,namex);
   Outtextxy(x1+10,y1-15,namey);
   Outtextxy(round((x1+x2)/2),y1-15,namegr);
   XMin:=X[1];
   XMax:=X[N];
   YMin:=1000;
   YMax:=-1000;
   for j:=1 to ngr do
   for i:= 1 to N do
      begin
         if(YMin>Y[j,i]) then YMin:=Y[j,i];
         if(YMax<Y[j,i]) then YMax:=Y[j,i];
      end;
   for i:=1 to ngr1 do
      begin
         Xtmp:= X1 + round((X2-X1)*(i-1)/(ngr1-1));
         Xt  := XMin + (XMax-XMin)*(i-1)/(ngr1-1);
         Str(Xt:6:2,XStr);
         OutTextXY(XTmp-35,y2+10,XStr);
         line(Xtmp,y1,Xtmp,y2);
      end;
   for i:=1 to ngr2 do
      begin
         Ytmp:=y1 + round((y2-y1)*(i-1)/(ngr2-1));
         Yt:=YMax - (YMax-YMin)*(i-1)/(ngr2-1);
         Str(Yt:6:2,YStr);
         OutTextXY(x1-55,YTmp,YStr);
         line(X1,Ytmp,X2,Ytmp);
      end;
   Setlinestyle(0,1,1);

   for j := 1 to ngr do
      begin
         SetColor(j);
         xs := x1 + trunc((x[1]-xmin)/(xmax-xmin)*(x2-x1));
         ys := y2 - trunc((y[j,1]-ymin)/(ymax-ymin)*(y2-y1));
         for i := 2 to n do
            begin
               xf := x1 + trunc((x[i]-xmin)/(xmax-xmin)*(x2-x1));
               yf := y2 - trunc((y[j,i]-ymin)/(ymax-ymin)*(y2-y1));
               line(xs,ys,xf,yf);
               xs:= xf;
               ys := yf;
            end;
      end;
   readln;
   CloseGraph;
end;
end.

Сообщение отредактировано: Гил Бейтс - 16.12.2007 21:53