Движения маятника Фуко по методу Эйлера

 

Program Mayatnik_Fuko;

const pi=3.1415; g=9.8; n=2000;

var m,gd,gm,I:integer;

r,w,t,tt,dt,v,d,L,vx0,vy0,x0,y0:real;

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

fuk:text;

begin

assign (fuk,'c:\fuko.DAT'); rewrite (fuk);

dt:=0.01; L:=50;

tt:=10;

t:=n*dt;

m:=round(tt/dt);

vx[0]:=0; vy[0]:=0;

x[0]:=1; y[0]:=0;

w:=0.04;

for i:=0 to m-1 do begin

R:=sqrt(x[I]*x[i]+y[i]*y[i]);

if r=0 then r:=1e-6;

v:=sqrt(vx[i]*vx[i]+vy[i]*vy[i]);

if v=0 then v:=1e-6;

vx[i+1]:=vx[i]+dt*(2*vy[i]*w+w*w*x[i]-g*x[i]/L);

vy[i+1]:=vy[i]+dt*(-2*vx[i]*w+w*w*y[i]-g*y[i]/L);

x[I+1]:=x[i]+dt*vx[i];

y[I+1]:=y[i]+dt*vy[i];

end;

writeln(y[m]);

for i:=0 to m do

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

close (fuk);

readln;

end.

Программа моделирования маятника Фуко

По исправленному методу Эйлера

 

Program Mayatnik_Fuko;

const pi=3.1415; g=9.8; n=2000;

var m,gd,gm,I:integer;

r,w,t,tt,dt,v,d,L,vxi1,vyi1,xi1,yi1:real;

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

fuk:text;

begin

assign (fuk,'c:\fukoisp1.DAT'); rewrite (fuk);

tt:=10;

dt:=0.01; L:=50;

m:=round(tt/dt);

t:=n*dt;

vx[0]:=0; vy[0]:=0;

x[0]:=1; y[0]:=0;

w:=0.04;

for i:=0 to m-1 do begin

R:=sqrt(x[I]*x[i]+y[i]*y[i]);

if r=0 then r:=1e-6;

v:=sqrt(vx[i]*vx[i]+vy[i]*vy[i]);

if v=0 then v:=1e-6;

vxi1:=vx[i]+dt*(2*vy[i]*w+w*w*x[i]-g*x[i]/L);

vyi1:=vy[i]+dt*(-2*vx[i]*w+w*w*y[i]-g*y[i]/L);

xi1:=x[i]+dt*vx[i];

yi1:=y[i]+dt*vy[i];

vx[i+1]:=vx[i]+0.5*dt*(2*vy[i]*w+w*w*x[i]-g*x[i]/L+2*vyi1*w+w*w*xi1-g*xi1/L);

vy[i+1]:=vy[i]+0.5*dt*(-2*vx[i]*w+w*w*y[i]-g*y[i]/L-2*vxi1*w+w*w*yi1-g*yi1/L);

x[I+1]:=x[i]+0.5*dt*(vx[i]+vxi1);

y[I+1]:=y[i]+0.5*dt*(vy[i]+vyi1);

end;

writeln(y[m]);

for i:=0 to m do

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

close (fuk);

readln;

end.

 

 

Колебания пружинного маятника

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

 

program avtokolebani;

const n=1500;

var I,d:integer;

h,w,b,k,v,v0,m,a,w1:real;

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

pm,vk:text;

begin

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

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

h:=0.01; a:=0.15; v0:=1000; k:=50; x[0]:=1; x[1]:=0;

w:=10;w1:=10;

b:=h*h*w*w;

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

begin

if I>100 then

c[i]:=h*k*sin((x[i]-x[I-1])/h)/(1+(x[i]-x[I-1])*(x[i]-x[I-1])/(v0*v0*h*h));

x[I+1]:=2*x[i]-x[I-1]-b+x[i]+c[i];

end;

for i:=1 to N do

begin

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

writeln (vk, i*h,' ',c[i]);

end;

close (PM); close (vk);

end.

Программа 2 по методу Эйлера

program avtokolebani;

const n=4000;

var I:integer;

h,w0,f,k,m,a,w1,d:real;

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

pm,vk:text;

begin

assign (pm,'c:\maitnik.DAT'); rewrite (pm);

assign (vk,'c:\maitnik1.DAT'); rewrite (vk);

h:=0.02; v[0]:=0;v[1]:=0;a:=0.005;

x[0]:=0.0; w0:=1.5;

w1:=1.5; f:=0.5;

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

begin

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

d:=-w0*w0*x[i]+f*cos(w1*i*h)-a*v[i]-w0*w0*x[i+1]+f*cos(w1*(i+1)*h)-a*v[i+1];

v[i+1]:=v[i]+0.5*h*d;

end;

for i:=1 to N do

begin

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

writeln (vk, i*h,' ',v[i]);

end;

close (PM); close (vk);

end.

Программа 3 по методу прогноза коррекции

 

program avtokolebani;

const n=2000;

var I,d:integer;

f,fi,alfa,delta,h,w,b,k,w0,m,a,w1:real;

v,x,v1,x1:array [0..n] of real;

pm,vk:text;

begin

assign (pm,'c:\mait05.DAT'); rewrite (pm);

assign (vk,'c:\diag05.DAT'); rewrite (vk);

h:=0.01;

a:=0.0;

k:=0.0;

f:=0.0;

alfa:=0;

delta:=0;

x[0]:=0.01; v[0]:=0;

w0:=3.1;

w1:=1.0;

x[0]:=1;v[0]:=0.0;w:=w0;

v[1]:=v[0]+h*(-w*w*x[0]+k*v[0]-a*v[0]-alfa*x[0]*x[0]*x[0]+delta*abs(x[0])*x[0]);

x[1]:=x[0]+h*v[0];

for i:=1 to N-1 do begin

v[i+1]:=v[i-1]+2*h*(-w*w*x[i]+k*v[i]-a*v[0]+f*cos(w1*i*h)-alfa*x[i]*x[i]*x[i]

+delta*abs(x[i])*x[i]);

x[i+1]:=x[i-1]+2*h*v[i];

fi:=-w*w*x[i+1]+k*v[i+1]-a*v[i+1]+f*cos(w1*(i+1)*h)-alfa*x[i+1]*x[i+1]*x[i+1]

+delta*abs(x[i+1])*x[i+1];

v1[i+1]:=v[i]+h/2*(-w*w*x[i]+k*v[i]-a*v[i]+f*cos(w1*i*h)-alfa*x[i]*x[i]*x[i]

+delta*abs(x[i])*x[i]+fi);

x1[i+1]:=x[i]+h/2*(v[i]+v[i+1]);

v[i+1]:=v1[i+1];

x[i+1]:=x1[i+1];

end;

for i:=1 to N do

begin

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

writeln (vk, x[i],' ',v[i]);

end;

close (PM); close (vk);

end.

Двойной маятник

 

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

program dm;

const n=5000;

var H,W,C,M,q:real;

i:integer; s:string;

x:array [1..N] of real;

y:array [1..N] of real;

f:text;

begin

assign (f,'d:\ft.dat');rewrite(f);

H:=0.1; W:=0.2; C:=2.5; Q:=0.02;

x[1]:=c;x[2]:=c;

y[1]:=0;y[2]:=0;

for i:=2 to n-1 do begin

x[i+1]:=2*x[i]-x[i-1]-H*H*w*w*(1+2*q)*x[i-1]+w*w*h*h*q*y[i-1];

y[i+1]:=2*y[i]-y[i-1]-H*H*w*w*y[i-1]+w*w*h*h*x[i-1];end;

for i:=1 to n do begin

writeln (f,i*h,x[i]+150);write (f,i*h,y[i]+100);end;

close (f);

end.

 

Связанные маятники

 

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

program ss;

const N=1000;

var h,w1,w2,k1,k2,a,B1,B2,C1,C2:real;

i:integer;

x1,x2:array[0..N] of real;

sm,sm1,sm2:text;

begin

w1:=1.1; w2:=1.09; h:=0.1; k1:=0.009; k2:=0.009; a:=0.0001;

assign(sm,'d:\maiat.dat');rewrite(sm);

assign(sm1,'d:\qmaiat.dat');rewrite(sm1);

assign(sm2,'d:\dmaiat.dat');rewrite(sm2);

B1:=2-h*h*w1*w1-h*h*k1-2*a*h;

B2:=2*h*a-1;

C1:=-h*h*w2*w2-h*h*k2+2;

C2:=h*h*k2;

x1[0]:=0.57;x1[1]:=0;

x2[0]:=0.0;x2[1]:=0;

for i:=1 to N-1 do begin

x1[i+1]:=x1[i]*B1+x1[i-1]*B2+x2[i]*k1;

x2[i+1]:=x2[i]*C1-x2[i-1]+x1[i]*C2;

end;

for i:=1 to N do

begin

writeln (sm,h*i,' ',x1[i]);

writeln (sm1,h*i,' ',x2[i]);

writeln (sm2,x1[i],' ',x2[i]);

end;

close (sm);

end.

Решение задачи Ферми-Паста-Улама

 

Здесь использован метод дискретного аналога ОДУ

 

program Fermi;

const nt=320; t=0.2; n=33;

var i,j:integer; k,k1,m,A,s,d:real;

y:array[1..nt,1..n] of real; E:array[150..150,1..n] of real;

mash,mash0,mash1:text;

begin

assign(mash0, 'c:\sash0.dat'); rewrite(mash0);assign(mash, 'c:\sash.dat');

rewrite(mash);assign(mash1, 'c:\sash2.dat'); rewrite(mash1);

k:=0.02; k1:=0.1; m:=100; A:=0;

for j:=1 to n do begin

y[1,j]:=0.1*sin(2*j*t); y[2,j]:=0.0;

end;

for i:=2 to nt-1 do begin

for j:=2 to n-1 do begin

s:=y[i-1,j+1]-y[i-1,j];

d:=y[i-1,j]-y[i-1,j-1];

y[i+1,j]:=t*t/m*(k*(y[i-1,j+1]-2*y[i-1,j])+k1*A*(s*s*s-d*d*d))+2*y[i,j]-y[i-1,j];

end;

writeln (i);

end;

for j:=2 to n-1 do begin

i:=300;

E[i,j]:=k*y[i,j]*y[i,j]/2+m/2*((y[i+1,j]-y[i,j])/t)*((y[i+1,j]-y[i,j])/t);

end;

for j:=1 to n do begin

writeln(mash0,j*t,' ',y[1,j]);

end;

i:=300;

for j:=2 to n-1 do begin

writeln(mash1,j*t,' ',E[i,j]);

end;

for j:=1 to n do begin

writeln(mash,j*t,' ',y[10,j],' ',y[20,j],' ',y[40,j],' ',y[60,j],' ',

y[80,j],' ',y[90,j],' ', y[100,j],' ',y[110,j],' ', y[120,j],' ',

y[130,j],' ', y[140,j],' ',y[150,j],' ', y[160,j],' ',y[170,j],' ',

y[180,j],' ',y[190,j],' ', y[200,j],' ',y[210,j],' ', y[220,j],' ',

y[230,j],' ', y[240,j],' ',y[250,j],' ', y[260,j],' ',y[270,j],' ',

y[280,j],' ',y[290,j],' ',y[300,j]);

end;

close(mash0); close(mash); close(mash1)

end.

Распространение волн на воде. Солитоны.