Движение тела
С учетом сопротивления среды
Движение небесного тела в гравитационном поле
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. Маятник Фуко
Программа для моделирования
|
