Скачиваний:
70
Добавлен:
03.10.2013
Размер:
35.84 Кб
Скачать

program dinam;

const m=2 (*количество решаемых ОДУ*); np=9; nk=6; nv=12;

type areal=array [1..m] of real;

var

del (* шаг*),k,x0 (*начальное значение аргумента *),x1 :real;

y0 (*начальные значения функций *),y1 :areal;

n (* число шагов *), kv (* кратность вывода *),i,kl,kl1,kl2 :integer;

vm,v,ak:array[1..nk] of real;

p:array[1..np] of real;

hg,h,s,pr:areal;

ro,pn,g,x:real;

label 10,20;

FUNCTION SGN(ab:real):integer;

begin

if ab>0 then sgn:=1;

if ab<0 then sgn:=-1;

if ab=0 then sgn:=0;

end;

procedure dydx(x:real;y:areal; var pr:areal);

(*расчет производных(pr) для текущих x,y*)

label 300,400;

var

vm:array[1..6] of real;

i:integer;

begin

g:=9.18;

h:=y;

p[8]:=pn * hg[1] / (hg[1] - h[1]);

p[9]:=pn*hg[2]/(hg[2]-h[2]);

p[6]:=p[8] + ro * g * h[1] *0.000001;

p[7]:=p[9]+ro*g*h[2]*0.000001;

v[1]:=ak[1] * SGN(p[1] - p[6]) * SQRt(ABS(p[1] - p[6]));

v[3]:=ak[3] * SGN(p[6] - p[3]) * SQRt(ABS(p[6] - p[3]));

v[5]:=ak[5] * sgn(p[7]-p[5])*sqrt(abs(p[7]-p[5]));

v[2]:=ak[2] * SGN(p[6] - p[2]) * SQRt(ABS(p[6] - p[2]));

v[4]:=ak[4] * SGN(p[7] - p[4]) * SQRt(ABS(p[7] - p[4]));

v[6]:=ak[6] * sgn(p[6] - p[7]) * sqrt(abs(p[6] - p[7]));

pr[1]:=(v[1]-v[2]-v[6]-V[3])/s[1];

pr[2]:=(v[6] - v[5] - v[4])/s[2];

FOR i:=1 TO 6 do

vm[i]:=v[i] * ro;

IF kl = 0 THEN goto 400;

IF kl = 1 THEN goTO 300;

300: writeln('x = ', x:nv:2, ' y[1]=',y[1]:nv:5,' y[2]=',y[2]:nv:5);

writeln('pr[1]=',pr[1],' pr[2]=',pr[2]);

400:

end;

procedure step(x:real; y:areal; del:real;

var x1:real;

var y1:areal);

(*шаг интегрирования:по X и Y расчет X1 и y1*)

var

pr,y12:areal;

i:integer ;

begin

dydx(x,y,pr);

for i:=1 to m do

y12[i]:=y[i]+del*pr[i]/2;

dydx(x+del/2,y12,pr);

x1:=x+del;

for i:=1 to m do

y1[i]:=y[i]+del*pr[i];

end;

begin

writeln('vvedite высотy емкостi hg[1]= /m/ ');

read( hg[1]);

writeln(' vvedite vicoty emkocti hg[2]= /m/ ');

read(hg[2]);

writeln('vvedite площад основаниa s[1]= /m^2/ ');

read(s[1]);

writeln('vvedite ploschad ocnovania s[2]= /m^2/ ');

read(s[2]);

writeln('vvedite плотность жидкости /кг/м^3/ = ');

read(ro);

writeln('vvedite начальное давление /МПА/ = ');

read(pn);

for i := 1 TO 5 do

begin

writeln('vvedite davlenie p[',i,']= /MPa/ ');

read(p[i]);

end;

for i := 1 TO 6 do

begin

writeln('vvedite koefficient cygaychego yctroictva ak[',i,']= ');

read(ak[i]);

end;

10: writeln('nachalnie yclovia');

writeln('vvedite x0 =');

read(x0);

writeln('vvedite y0[1] = ');

read(y0[1]);

writeln('vvedite y0[2] = ');

read(y0[2]);

20: writeln('vvedite число шагов n =');

read(n);

цкшеудт(эммувшеу краттость вывода лм =э)ж

read(kv);

writeln('vvedite del');

read(del);

writeln('вывод промежуточ.: нет - 0 частич. - 1 полн.- 2 ');

read(kl1);

writeln('вывод результатов: частич. - 1 полн.- 2 ');

read(kl2);

for i:=1 to n do

begin

kl:=kl1;

step(x0,y0,del,x1,y1);

x0:=x1;

y0:=y1;

if (i div kv) =(i/kv) then

begin

write('ПРИ i=',i:3);

if kl2=1 then

writeln(' x=',x0:10:5,' y[1]=',y0[1]:nv:5,' y[2]=',y0[2]:nv:5);

if kl2=2 then

begin

kl:=kl2;

dydx(x0,y0,pr);

end

end

end;

writeln('прекращение счета - 0; продолжение - 1; счет с новыми н.у.- 2 ')

read(kl);

if kl=1 then

goto 20;

if kl=2 then

goto 10;

end.

Соседние файлы в папке Программы (Pascal) - 2003