A Maple egy fejlett matematikai problémamegoldó és programozói szoftver.
-
xDJCx #441
A Kepleres feladathoz kapcsolódva az alábbi procedura egy bolygó Nap körüli mozgását szimulálja Newton-törvény és a gravitációs törvény alapján.
Az elméleti háttere Feynman: Mai fizika I. kötetében le van írva, és onnan származnak a példa kezdeti paraméterei értékei is. (Lehetne szebben is megírni, de most annyi időt nem töltöttem vele.)
restart;with(plots):
bolygo:=proc(N,dt,x0,y0,vx0,vy0,M,G)
description "Bolygómozgás számítása Gravitációs és Newton törvénybõl Forrás: Feynmann : Mai fizika I. ";
#N;# iterációk száma
#dt;# idõlépésköz
#x0;y0;vx0;vy0;# t=0-beli kezdeti pálya adatok
#M,G;# Nap tömege és gravitációs állandó
local x,y,xyuj,vx,vy,ax,ay,vxtmdt2,vytmdt2,vxtpdt2,vytpdt2,txy,nap,k,r,Txy,xt,yt,axt,ayt,bolygo,mozgo,ujpalyaadatok;
use plots in
x:=Array(1..N):y:=Array(1..N):vx:=Array(1..N):vy:=Array(1..N):ax:=Array(1..N):ay:=Array(1..N):
txy:=Array(1..N):
xyuj:=proc(xt,yt,vxtmdt2,vytmdt2)
local r,axt,ayt,xuj,yuj,vxtpdt2,vytpdt2:
r:=sqrt(xt^2+yt^2):
axt:=-G*M*xt/r^3:ayt:=-G*M*yt/r^3:
vxtpdt2:=vxtmdt2+dt*axt:
vytpdt2:=vytmdt2+dt*ayt:
xuj:=xt+dt*vxtpdt2:
yuj:=yt+dt*vytpdt2:
[xuj,yuj,vxtpdt2,vytpdt2]:
end proc;
txy[1]:=[0,x0,y0]:#t=0
xt:=x0:yt:=y0:
r:=sqrt(xt^2+yt^2):# t=0-ban
axt:=-G*M*xt/r^3:ayt:=-G*M*yt/r^3:# t=0-ban
vxtmdt2:=vx0+dt/2*axt:vytmdt2:=vy0+dt/2*ayt: # t=dt/2-ban
# t=dt
xt:=x0+vxtmdt2*dt;yt:=y0+vytmdt2*dt:#t=dt
txy[2]:=[dt,xt,yt]:# t=dt
ujpalyaadatok[1]:=xt:ujpalyaadatok[2]:=yt:
ujpalyaadatok[3]:=vxtmdt2:ujpalyaadatok[4]:=vytmdt2:
for k from 3 to N do
xt:=ujpalyaadatok[1];yt:=ujpalyaadatok[2];
vxtmdt2:=ujpalyaadatok[3]:vytmdt2:=ujpalyaadatok[4]:
ujpalyaadatok:=xyuj(xt,yt,vxtmdt2,vytmdt2):
txy[k]:=[(k-1)*dt,ujpalyaadatok[1],ujpalyaadatok[2]]:
end do:
Txy:=seq([txy[t][2],txy[t][3]], t=1..N):
bolygo:= proc(xy) plots[pointplot]([xy],color=blue,symbol=solidcircle,symbolsize=20) end proc:
mozgo:=seq(bolygo(Txy[n]), n=1..N):nap:=plots[pointplot]([0,0],color=yellow,symbol=solidcircle,symbolsize=40):
display(nap,display(mozgo,scaling=constrained,insequence = true,axes=normal));
end use;
end proc;
bolygo(82,0.05,0.5,0,0,1.63,1,1);# G és M 1-nek választva itt.