Я записал исходные данные координат точек образования сопла, построил его, и нашел число маха 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(c));
Yb[2,i] := exp(-1.0/(gamma-1)*ln(c));
Yb[3,i] := exp(-1.0/2.0*ln(c));
Yb[4,i] := exp(-1.0*ln(c));
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(c));
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('ў ¤ (r)¬ ¤Ё Ї §(r)Ґ Є(r)аҐм Ґбвм');}
{
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('ў ¤ (r)¬ ¤Ё Ї §(r)Ґ Є(r)ап Ґв');
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.