# Newton's and Kepler's laws# Equlibrium fo Two Forcesrestart;R1 := G*m[1]*m[x]/rx^2 = G*m[x]*m[2]/(r - rx)^2;Sr := solve(R1, rx);rx := factor(Sr[2]):rnx := simplify(subs(r = 1, m[2] = N^2*m[1], rx),symbolic):
rnx := subs(N = sqrt(n), rnx);plot(rnx, n=0..10, labels=[`n`, 'rx']);# Equlibrium fo Three Forcesr1 := sqrt((x1 - x)^2 + y^2):
r2 := sqrt((x2 - x)^2 + y^2):
r3 := sqrt(x^2 + (y3 - y)^2):
F1 := G*m1*m/r1^2:
F2 := G*m2*m/r2^2:
F3 := G*m3*m/r3^2:
F1x := F1*(x1 - x)/r1: F1y := -F1*y/r1:
F2x := F2*(x2 - x)/r2: F2y:= -F2*y/r2:
F3x := -F3*x/r3: F3y := F3*(y3 - y)/r3:Rx := F1x + F2x + F3x=0; Ry := F1y + F2y + F3y=0;System := subs(G = 1, m = 1, m1 = 12, x1 = -5, m2 = 9,
x2 = 3, m3 = 6, y3 = 8, {Rx, Ry});with(plots):
implicitplot(System, x = -4..2, y = 0..6, grid=[60,60],
labels = [`x`,`y`], color = black);Digits:=20:
NS1 := fsolve(System, {x, y}, {x = -1..0, y = 0..1});
NS2 := fsolve(System, {x, y}, {x = -1..0, y = 4..5});evalf(subs(NS1, System)); evalf(subs(NS2, System));# Equlibrium fo Three Forces Computed from the Potential EnergyDigits:=10:
U := - G*m1/r1 - G*m2/r2 - G*m3/r3:
with(linalg):
g := grad(U,[x,y]);Un := subs(G = 1, m = 1, m1 = 12, x1 = -5,
m2 = 9, x2 = 3, m3 = 6, y3 = 8, U);gn := grad(Un, [x, y]);implicitplot({seq(Un = -i/6, i = 24..50)}, x=-6..3, y=-1..9,
scaling = constrained, labels = [`x`, `y`], color = black,
grid=[60,60]);
plot3d(Un, x = -6..3, y = -1..9, view = -12..-2,
orientation = [-100, 75], style = hidden, color = black,
grid=[60,60], axes = boxed, labels = [`x`, `y`, 'U']);NS1 := fsolve({gn[1], gn[2]}, {x, y}, {x = -1..0, y= 0..1});
NS2 := fsolve({gn[1], gn[2]}, {x, y}, {x = -1..0, y= 4..5});subs(NS1, hessian(Un, [x, y]));
subs(NS2, hessian(Un, [x, y]));# Gravitation of the Massive Line Segmentrestart;
with(plots): with(linalg):
r := sqrt((xi - x)^2 +y^2 +z^2):
sigma := m/L: # linear mass density
U := Int(G*sigma/r, xi = -L/2..L/2);U := factor(value(U));Us := subs(G=1,m=1,L=1, z=1/5, U);plot3d(Us, x = -2..2, y=-2..2, axes = boxed, style = hidden,
color = black, orientation = [-80, -130],
grid=[35,35], labels = [`x`, `y`, 'U']);
p1 := implicitplot({seq(Us = i/10, i = 5..30)}, x = -1..1,
y = -1..1, color = black, scaling = constrained,grid=[40,40]):p2 := gradplot(Us, x = -1..1, y = -1..1, color = black,
scaling = constrained, arrows = SLIM, grid=[15,15]):
display({p1,p2},labels = [`x`,`y`]); Uinf := Limit(U, L=0);Uinf := value(Uinf);ginf := grad(Uinf, [x, y, z]);abs(ginf) = simplify(norm(ginf, 2), symbolic);Us := subs(G = 1, m = 1, L = 1, U):
D2r := [diff(x(t), t, t), diff(y(t), t, t), diff(z(t), t, t)]:
g := subs(x = x(t), y = y(t), z = z(t), grad(Us, [x, y, z])):
IniC:= x(0) = 1, D(x)(0) = 0, y(0) = 0, D(y)(0) = 1,
z(0) = 3/4, D(z)(0) = 0:
Ns := dsolve({seq(D2r[i] = g[i], i = 1..3), IniC},
{x(t), y(t), z(t)}, numeric);for i from 0 to 500 do;
T := i/12.5;
NsT := Ns(T):
X[i] := subs(NsT,x(t)); Vx[i] := subs(NsT,diff(x(t),t));
Y[i] := subs(NsT,y(t)); Vy[i] := subs(NsT,diff(y(t),t));
Z[i] := subs(NsT,z(t)); Vz[i] := subs(NsT,diff(z(t),t));
KepVec[i] := convert(crossprod([X[i], Y[i] ,Z[i]],
[Vx[i], Vy[i], Vz[i]]), list)[2..3];
KepAbs[i] := norm(KepVec[i], 2);
od:spacecurve({[seq([X[i], Y[i], Z[i]], i = 0..500)],
[[-1/2, 0, 0], [1/2, 0, 0]]}, labels=[`x`, `y`, `z`],color=black,axes=boxed);plot([seq(KepVec[i], i = 0..500)], labels=[`y`, `z`],color=black);plot([seq([i/25, KepAbs[i]], i = 0..500)],
labels = ['t', 'MofI']);IniC := x(0) = 1, D(x)(0) = 0, y(0) = 0, D(y)(0) = 1,
z(0) = 1/2, D(z)(0) = 1/2:# The Earth Satelliterestart;
dx := diff(x(t), t, t) = -G*Mz*x(t)/(x(t)^2 + y(t)^2)^(3/2):
dy := diff(y(t), t, t) = -G*Mz*y(t)/(x(t)^2 + y(t)^2)^(3/2):
G:=6.67*10^(-11): Mz:=6*10^24:
Inic := x(0) = 7*10^(6), D(x)(0)=0, y(0)=0, D(y)(0)=9*10^3:
Digits := 15:
Ns := dsolve({dx, dy, Inic}, {x(t),y(t)}, numeric):for i from 0 to 400 do;
T := i*40;
NsT := Ns(T):
X[i] := subs(NsT,x(t)); Vx[i] := subs(NsT,diff(x(t),t));
Y[i] := subs(NsT,y(t)); Vy[i] := subs(NsT,diff(y(t),t));
MofI[i] := X[i]*Vy[i]-Y[i]*Vx[i];
od:
with(plots):
p1 := polarplot(6378*10^3, phi = 0..2*Pi):
p2 := plot([seq([X[i], Y[i]], i = 0..327)], thickness=2):
display({p1, p2}, labels = ['x', 'y'], scaling = constrained);
plot([seq([i*40,MofI[i] - 0.63*10^11], i = 0..400)],
labels=['t','Wp']);tx := evalf((t1+t2)/2):
t1 := 326*40: t2 := 327*40:
y1 := subs(Ns(t1),y(t));
y2 := subs(Ns(t2),y(t));while (t1 < tx) and (tx < t2) do:
yx := subs(Ns(tx),y(t)):
if yx > 0 then
y2 := yx: t2 := tx:
else
y1 := yx: t1 := tx:
fi;
od:
Tx = evalf(tx);Hx := floor(tx/3600):
Mx := floor((tx - Hx*3600)/60):
Sx := tx - Hx*3600 - Mx*60:
Hx, Mx, Sx;XS := [seq(X[i], i = 0..328)]: YS := [seq(Y[i], i = 0..328)]:
VxS := [seq(Vx[i], i = 0..328)]: VyS := [seq(Vy[i], i = 0..328)]:
save(G, Mz, XS, YS, VxS, VyS,`orbit.sav`);# The Earth Satellite, Second Solutionrestart;with(plots):read(`orbit.sav`):ax := -G*Mz*x/(x^2 + y^2)^(3/2):
ay := -G*Mz*y/(x^2 + y^2)^(3/2):
i := 'i': j := i + 1:
for k from 0 to 3 do:
x := 7*10^6: Vx := 0:
y := 0: Vy := 9000:
dt := evalf(1/2^k);
for i from 0 to 328 do:
X[i] := evalf(x); Y[i] := evalf(y):
for n from 1 to 40*2^k do:
x := evalf(ax*dt^2/2 + Vx*dt + x);
y := evalf(ay*dt^2/2 + Vy*dt + y);
Vx := evalf(ax*dt + Vx);
Vy := evalf(ay*dt + Vy);
od:
if i mod 41 = 0 then
dX[k,i]:=X[i]-XS[j]; dY[k,i]:=Y[i]-YS[j];
fi;
od: p[k] := plot([seq([(X[i] - XS[j])/1000,
(Y[i] - YS[j])/1000], i = 0..328)], color = black):
od:p1 := display({seq(p[k], k = 0..3)}, thickness = 3):p1;SI := [seq(i*41, i = 0..8)]:
p2 := plot({seq([seq([dX[k, i]/1000, dY[k, i]/1000],
k = 0..1), [0, 0]], i = SI)}, color = black):
display({p1, p2}, scaling = constrained, labels = ['dx', 'dy']);
display({p1, p2}, view = [-0.1..0.5, -0.4..0.2],
scaling = constrained, labels = ['dx', 'dy']);# The Lost Screwrestart;read(`orbit.sav`):
dx := diff(x(t), t, t) = -G*Mz*x(t)/(x(t)^2 + y(t)^2)^(3/2):
dy := diff(y(t), t, t) = -G*Mz*y(t)/(x(t)^2 + y(t)^2)^(3/2):
Inic := x(0) = 7*10^(6), D(x)(0)=1, y(0)=0, D(y)(0)=9*10^3:
Ns := dsolve({dx, dy, Inic}, {x(t),y(t)}, numeric): for i from 0 to 328 do;
T := i*40; NsT := Ns(T);
X[i] := subs(NsT,x(t)); Vx[i] := subs(NsT,diff(x(t),t));
Y[i] := subs(NsT,y(t)); Vy[i] := subs(NsT,diff(y(t),t));
od:i := 'i': j := i + 1:
plot([seq([X[i] - XS[j], Y[i] - YS[j]], i = 0..328)],
labels = [`d x`, `d y`]);
plot([seq([Vx[i] - VxS[j], Vy[i] - VyS[j]], i = 0..328)],
labels = [`d Vx`, `d Vy`]); R := sqrt(XS[j]^2 + YS[j]^2):
Sin := YS[j]/R: Cos:=X[j]/R:
GCX := (X[i] - XS[j])*Cos - (Y[i] - YS[j])*Sin:
GCY := (X[i] - XS[j])*Sin + (Y[i] - YS[j])*Cos:
GCVx := (Vx[i] - VxS[j])*Cos - (Vy[i] - VyS[j])*Sin:
GCVy := (Vx[i] - VxS[j])*Sin + (Vy[i] - VyS[j])*Cos:plot([seq([GCX, GCY], i = 0..328)], labels = [`d x`, `d y`]);
plot([seq([GCVx, GCVy], i = 0..327)], labels = [`d Vx`,`d Vy`]);