Движение тела

С учетом сопротивления среды

 

Program a; Const n=5000;g=9.8;M=60;h=0.2; Label 1; Var c,S,q,k2:real; i,p,r:integer; X,V:array [0..n] of real; Dnt:text; Begin Assign(dnt,'c:\ayrat.dat'); Rewrite(dnt); x[0]:=n*h;v[0]:=0;c:=0.44;S:=0.9; q:=1.29; k2:=0.5*c*S*q; for i:=0 to n-1 do

Begin v[i+1]:=h*g-k2*h*v[i]*v[i]/m+v[i]; x[i+1]:=-h*v[i]+x[i]; if x[i]<200 then begin s:=20;c:=0.4;q:=1.;k2:=0.5*c*s*q;end; if x[i]<0 then begin x[i]:=0;v[i]:=0;end; if x[i]=0 then begin x[i+1]:=0;v[i+1]:=0;end; end; For i:=0 to n do begin writeln(dnt,i*h,' ',x[i],' ',v[i]); {writeln(dnt,i,' ',R[i]);} end; close(dnt); end.

Движение небесного тела в гравитационном поле

 

program grav;

const n=2000;G=6.67;M=2;h=0.033;

var i:integer;

f1,f2,v11,v21,x1,y1,a,v0,alfa:real;

X,Y,V1,V2,R:array [0..n] of real;

dnt:text;

Begin

assign(dnt,'c:\ast1.dat');

rewrite(dnt);

a:=30;

alfa:=a*3.14/180; v0:=0.7;

x[0]:=-10;y[0]:=0;v1[0]:=v0*cos(alfa);v2[0]:=v0*sin(alfa);

for i:=0 to n-1 do

begin

f1:=G*M*x[i]/(sqrt(x[i]*x[i]+y[i]*y[i])*(x[i]*x[i]+y[i]*y[i]));

v11:=v1[i]-h*f1;

f2:=G*M*y[i]/(sqrt(x[i]*x[i]+y[i]*y[i])*(x[i]*x[i]+y[i]*y[i]));

v21:=v2[i]-h*f2;

x1:=x[i]+h*v11;

y1:=y[i]+h*v21;

v1[i+1]:=v1[i]-0.5*h*(f1+G*M*x1/(sqrt(x1*x1+y1*y1)*(x1*x1+y1*y1)));

v2[i+1]:=v2[i]-0.5*h*(f2+G*M*y1/(sqrt(x1*x1+y1*y1)*(x1*x1+y1*y1)));

x[i+1]:=x[i]+0.5*h*(v11+v1[i]);

y[i+1]:=y[i]+0.5*h*(v21+v2[i]);

end;

For i:=0 to n do

begin

writeln(dnt,x[i],' ',y[i]);

end;

close(dnt);

end.

 

Полет сверхзвукового самолета

 

program Polet;

const N=1000;

var i:integer;

x,y,v,a,vx,vy1:array[0..N] of real;

xy,vy,ay:TEXT;

h,k1,k2,r,s,k,T,m,g,r0,v0,a0,a1,Ft,Nr,Q:real;

begin

assign(xy,'C:\polet_xy1.dat');rewrite(xy);

assign(vy,'C:\polet_vy2.dat');rewrite(vy);

assign(ay,'C:\polet_ay3.dat');rewrite(ay);

g:=9.8; m:=2000; v0:=100; a0:=45;

a1:=3.1415*a0/180;

h:=0.05; k1:=0.02; k2:=0.2; r0:=0.01;

k:=138; s:=50; T:=300; v[0]:=v0;

a[0]:=a1;

y[0]:=0;

for i:=0 to N-1 do

begin

vx[i]:=v[i]*cos(a[i]);vy1[i]:=v[i]*sin(a[i]);

r:=r0{*exp(-m*g*(i*h)/(k*T))};

Ft:=(1+0.5*exp(-sqr(0.001*v[i]-1.5)))*(10000-0.27*y[i]);

Q:=0.5*k1*r*s*sqr(v[i]);

Nr:=0.5*k2*r*s*sqr(v[i]);

x[i+1]:=x[i]+vx[i]*h;

y[i+1]:=y[i]+vy1[i]*h;

v[i+1]:=v[i]+h*(Ft-Q-m*g*cos(a[i]))/(m*v[i]*sin(a[i]));

a[i+1]:=a[i]+h*(Nr-m*g*cos(a[i]))/(m*sqr(v[i])*sin(a[i]));

end;

for i:=0 to N do

begin

writeln(xy,x[i],' ',y[i]);

writeln(vy,y[i],' ',v[i]);

writeln(ay,y[i],' ',a[i]);

end;

close(xy);

close(vy);

close(ay);

end.

 

Полет одноступенчатой и баллистической ракеты

program Ballistic_raket;

var i,n:integer;

x,y,v,h,k1,k2,alpha,m0,m,mk,b,g0,g,c,s:real;

R,Ln10,u,Ft,Fs,ro0,teta,Vzvuk,a,tet:real;

FileXYM:text;

 

Procedure NextXi;

begin

g:=g0*R*R/((R+y)*(R+y));

Ft:=u*alpha;

Fs:=k1*v+k2*exp(-b*y)*ln10*v*v;

Fs:=k2*exp(-b*y)*ln10*v*v;

a:=(Ft-Fs)/m;

if v<=Vzvuk then begin y:=y+h*v;v:=v+h*(a-g);

end

else

begin {SinUg:=sin(3.14/2+teta);

CosUg:=cos(3.14/2+teta);}

y:=y+h*V*Sin(teta);x:=x+h*v*Cos(teta);

Teta:=teta-(h*g*Cos(teta)/v);

v:=v+h*(a-g*Sin(teta));

end;

writeln(FileXYM,x,' ',y,' ',v,' ',a,' ',m,' ',h*i);

end;

begin

assign(FileXYM,'ballist2.dat');rewrite(FileXYM);

n:=5000;b:=0.0000125;

ln10:=ln(10);s:=0.00005;c:=2;k2:=s*c*0.5;

m0:=20;ro0:=1300000;mk:=5;

alpha:=0.3;R:=6400;u:=1.05;h:=0.1;

g0:=0.0098;

VZvuk:=0.34;

v:=0;

x:=0;y:=0;

tet:=40;

teta:=tet*3.1415/180;

m:=m0;i:=0;

while (y>=0) and (i<5000) do

begin

if m>mk then m:=m0-alpha*h*i;

if m<mk then m:=mk;

NextXi;

Inc(i);

end;

Close(FileXYM);

end.

Стыковка космического корабля

program stykovka

Uses Graph,crt;

const n=1000;

Var grDriver, grMode, ErrCode: integer;

ch:char;

x,y,vx,vy,dx,dy:array [0..n] of real;

i:integer;

h,w:real;

dnt:text;

Begin

assign(dnt,'d:graf12.dat');

rewrite(dnt);

grDriver:=Detect;

InitGraph(grDriver, grMode, '');

ErrCode:=GraphResult;

If ErrCode <> grOk Then Halt(1);

setbkcolor(1);

x[0]:=10;y[0]:=10;vx[0]:=0;vy[0]:=0;dx[0]:=0;dy[0]:=0;

h:=0.33;w:=0.02;i:=0;

begin

setcolor(4);

setfillstyle(solidfill,2);

bar(400,300,440,320);

repeat

repeat

for i:=0 to n-1 do

begin

if keypressed then ch:=readkey;

case ch of

#72:begin dy[i]:=dy[i]-0.5;

end;

#75:begin dx[i]:=dx[i]-0.5;

end;

#77:begin dx[i]:=dx[i]+0.5;

end;

#80:begin dy[i]:=dy[i]+0.5;

end;

#27: halt;

end;

begin

x[i+1]:=x[i]+h*vx[i];

y[i+1]:=y[i]+h*vy[i];

vx[i+1]:=vx[i]+h*(3*w*w*y[i]-2*vy[i]*w+dx[i]);

vy[i+1]:=vy[i]+h*(2*vx[i]*w+dy[i]);

end;

repeat

begin

setcolor(7);

circle(round(x[i]),round(y[i]),5);

begin

writeln (dnt,x[i],'',y[i]);

end;

delay(15000);

setcolor(1);

circle(round(x[i]),round(y[i]),5);

end;

until keypressed;

end;

until ((x[i]>400) and (x[i]<440));

until ((y[i]>300) and (y[i]<320));

for i:=0 to n do

begin

setcolor(13);

OutTextXY(10,10,' стыковка');

readln;

end;

closegraph;

close(dnt);

end;

end.

6.7. Параметрический маятник

Программа 1 по методу дискретного аналога

 

program PM1;

const g=9.81;l=5;n=1500;

var i:integer;

h,w0,w1,t0,t1,m,a,l1:real;

f,x:array [0..n] of real;

y:array [0..n] of real;

PARM:text;

begin

h:=0.1; a:=0.5; t0:=1; t1:=20; m:=30;l1:=1;

x[0]:=1; x[1]:=0; y[0]:=l*cos(a*3.14/180); y[1]:=0;

f[0]:=0.01; f[1]:=0.01;

w0:=sqrt(g/l);w1:=2*w0;

assign (PARM,'d:\maitnik.DAT'); rewrite (PARM);

for i:=1 to (N-1) do

begin

f[i+1]:=2*f[i]-f[i-1]-h*h*g*sin(f[i])/(l+l1*cos(w1*i*h));

x[i]:=sin(f[i])*(l+l1*cos(w1*i*h));

y[i]:=cos(f[i])*(l+l1*cos(w1*i*h));

end;

for i:=1 to N do

begin

writeln (PARM,i*h,' ',x[i],' ',y[i]);

end;

close (PARM);

end.

 

Программа 2 с понижением порядка ОДУ

и использованием исправленного метода Эйлера

 

program PM1;

const g=9.81;l=5;n=1000;

var i:integer;

f1,u1,h,w0,w1,t0,t1,m,a,l1:real;

f,x,u:array [0..n] of real;

y: array [0..n] of real;

PARM:text;

begin

h:=0.05;l1:=1;

x[0]:=0;y[0]:=0;f[0]:=0.5;u[0]:=0.0;

w0:=sqrt(g/l);w1:=2*w0;

assign(PARM,'c:\m1.DAT');rewrite(PARM);

for i:=0 to (N-1) do

begin

f1:=f[i]+h*u[i];

u1:=u[i]-h*g*sin(f[i])/(l+l1*cos(w1*i*h));

f[i+1]:=f[i]+0.5*h*(u[i]+u1);

u[i+1]:=u[i]-0.5*h*g*(sin(f[i])/(l+l1*cos(w1*i*h))

+sin(f1)/(l+l1*cos(w1*h*(i+1))));

x[i]:=sin(f[i])*(l+l1*cos(w1*i*h));

y[i]:= cos(f[i])*(l+l1*cos(w1*i*h));

end;

for i:=0 to N do

begin

writeln (PARM,i*h,' ',x[i]);

end;

close(PARM);

end.

Маятник Фуко

 

Программа для моделирования