Harjoitus 11 -- Ratkaisut

1

Tässä harjoituksessa ei käsitellä kahden pisteen reuna-arvoprobleemaa. Sitä käsitellään 3. harjoitustyössä ja harjoitus 11:ssa. Tässä harjoituksessa lähinnä esitellään differentiaaliyhtälöiden ratkaisumahdollisuuksia Mathematicalla.

Maan säde, massa ja gravitaatiovakio löytyvät paketista PhysicalConstants. Supistetaan yksiköt pois.

<<PhysicalConstants`

R=EarthRadius/Meter
M=EarthMass/Kilogram
G=GravitationalConstant/(Meter^2)/Newton*(Kilogram^2)

Differentiaaliyhtälöt paikkavektorin r komponenteille:

eqx=x''[t]==-((G M)/(((((x[t])^2)+((y[t])^2)))^(3/2)))x[t]-(γ/m)(Sqrt[((x'[t])^2)+((y'[t])^2)])x'[t];
eqy=y''[t]==-((G M)/(((((x[t])^2)+((y[t])^2)))^(3/2)))y[t]-(γ/m)(Sqrt[((x'[t])^2)+((y'[t])^2)])y'[t];

a) ei kitkavoimaa

Ilmaistaan lähtöpaikka leveysasteen avulla:

θ=60.0Degree; (*lähtöpaikan leveysaste*)
r0=R{Cos[θ],Sin[θ]};

Lähtöpaikka voidaan piirtää ympyrän kehälle:

Show[Graphics[Circle[{0,0},R]],Graphics[Point[r0]],Axes->True]

Olkoon α lähtönopeuden ja Maan pinnan välinen kulma lähtöpaikassa. Tällöin nopeuden suuntakulma sellaisessa koordinaatistossa, jossa y-akseli osoittaa Maan akselin suuntaisesti, on ϕ=π/2-(θ+α) (osoita!)

α=20.0Degree;
V=10000; (*lähtönopeus metreinä sekunnissa*)
ϕ=(θ+α)-π/2;
v0=V{Cos[ϕ],Sin[ϕ]};

Vektoria (Subscript[v, 0]) on hieman suurenneltava, jotta se näkyisi Maan mittakaavassa.

g0=Show[Graphics[Circle[{0,0},R]],Graphics[Point[r0]],Graphics[Arrow[{r0,r0+200v0}]],Axes->True]

Lähdetään ratkaisemaan yhtälöitä. Kannattaa aloittaa pienestä T:n arvosta (jos esim. kappaleen lentorata vie Maan keskipisteeseen, NDSolve kaatuu).

T=18000;(*lentoaika sekunteina*)
val={m->10000,γ->0.0};
sol=Flatten[NDSolve[{eqx,eqy,x[0]==r0[[1]],y[0]==r0[[2]],x'[0]==v0[[1]],y'[0]==v0[[2]]}/.val,{x[t],y[t]},{t,0,T}]];
{X[t_],Y[t_]}={x[t],y[t]}/.sol;
Show[g0,ParametricPlot[{X[t],Y[t]},{t,0,T},AspectRatio->1]]

Ilman vastusvoimia kappaleen lentorata on ellipsi, jonka toinen polttopiste on Maan keskipisteessä. Jotta kappale palaisi lähtöpisteeseensä, on lähtönopeuden oltava Maan pinnan suuntainen (muuten kappale törmäisi Maan pintaan jossain muualla kuin lähtöpisteessä -- laskennallisesti kappale voi tietenkin matkata alueessa jossa etäisyys Maan keskipisteestä on <R). Kirjoitetaan ratkaisuprosessi omaksi funktiokseen, jonka argumenit ovat
	θ	lähtöpaikan leveysaste
	α	lähtönopeuden suuntakulma
	V	lähtönopeus
	m	kappaleen massa
	γ	kitkakerroin
	T	lentoaika

raketti[θ_,α_,V_,m0_,γ0_,T_]:=
Module[{ϕ,r0,v0,X,Y,eqx,eqy,g0},
ϕ=θ+α-π/2;
r0=R{Cos[θ],Sin[θ]};
v0=V{Cos[ϕ],Sin[ϕ]};
g0=Show[Graphics[Circle[{0,0},R]],Graphics[Point[r0]],Graphics[Arrow[{r0,r0+200v0}]],Axes->True];
eqx=x''[t]==-((G M)/(((((x[t])^2)+((y[t])^2)))^(3/2)))x[t]-(γ/m)(Sqrt[((x'[t])^2)+((y'[t])^2)])x'[t];
eqy=y''[t]==-((G M)/(((((x[t])^2)+((y[t])^2)))^(3/2)))y[t]-(γ/m)(Sqrt[((x'[t])^2)+((y'[t])^2)])y'[t];
val={m->m0,γ->γ0};
sol=Flatten[NDSolve[{eqx,eqy,x[0]==r0[[1]],y[0]==r0[[2]],x'[0]==v0[[1]],y'[0]==v0[[2]]}/.val,{x[t],y[t]},{t,0,T}]];
{X[t_],Y[t_]}={x[t],y[t]}/.sol;
Show[g0,ParametricPlot[{X[t],Y[t]},{t,0,T},AspectRatio->1]]
];

raketti[60.0Degree,0.0Degree,10000,10000,0.0,20000]

b) joo kitkavoimaa

Otetaan käyttöön epärealistinen malli kitkavoimalle (kitkakerroin on sama riippumatta etäisyydestä Maan pinnasta).

raketti[60.0Degree,20.0Degree,15000,15000,0.001,6500]

Raketti ei plaaa lähtöpaikalleen millään parametrien arvoilla (kitkavoima kääntää sen kulkusuunnan aina Maan keskipistettä kohti).

b-ekstra) parempi kitkavoima

Muutetaan hieman kitkavoiman lauseketta, muodostetaan voima joka heikkenee kun raketti loittonee Maasta. Heikkeneminen voisi olla vaikka muotoa 1/(r(t)/R) (täysin hatusta vedetty).

raketti2[θ_,α_,V_,m0_,γ0_,T_]:=
Module[{ϕ,r0,v0,X,Y,eqx,eqy,g0},
ϕ=θ+α-π/2;
r0=R{Cos[θ],Sin[θ]};
v0=V{Cos[ϕ],Sin[ϕ]};
g0=Show[Graphics[Circle[{0,0},R]],Graphics[Point[r0]],Graphics[Arrow[{r0,r0+200v0}]],Axes->True];
eqx=x''[t]==-((G M)/(((((x[t])^2)+((y[t])^2)))^(3/2)))x[t]-(γ/m)((Sqrt[((x'[t])^2)+((y'[t])^2)])/((Sqrt[((x[t])^2)+((y[t])^2)])/R))x'[t];
eqy=y''[t]==-((G M)/(((((x[t])^2)+((y[t])^2)))^(3/2)))y[t]-(γ/m)((Sqrt[((x'[t])^2)+((y'[t])^2)])/((Sqrt[((x[t])^2)+((y[t])^2)])/R))y'[t];
val={m->m0,γ->γ0};
sol=Flatten[NDSolve[{eqx,eqy,x[0]==r0[[1]],y[0]==r0[[2]],x'[0]==v0[[1]],y'[0]==v0[[2]]}/.val,{x[t],y[t]},{t,0,T}]];
{X[t_],Y[t_]}={x[t],y[t]}/.sol;
Show[g0,ParametricPlot[{X[t],Y[t]},{t,0,T},AspectRatio->1]]
];

raketti2[60.0Degree,20.0Degree,15000,15000,0.001,13000]

raketti2[60.0Degree,0.0Degree,15000,12000,0.001,4000]

b-ekstra2)

Lisäksi voitaisiin lisätä yhtälöihin työntövoima, esim.
	f(t)={({{(Subscript[f, 0])v(t), kun polttoainetta on jäljellä}, {0, kun polttoaine on loppu}})
Nyt ratkaisussa olisi määrättävä polttoaineen määrä alussa (Subscript[m, t]) ja polttoaineen palonopeus c (kg/s). Kappaleen massa olisi myös ajan funktio:
	m(t)={({{(Subscript[m, 0])-c t, kun t<(Subscript[m, t])/c}, {(Subscript[m, 0])-(Subscript[m, t]), kun t≥(Subscript[m, t])/c}})
(Subscript[m, t])/c on siis se ajanhetki, kun polttoaine loppuu.

2

Muodostetaan differenssiyhtälöt. Ilman kitkatermiä niistä saadaan eksplisiittisesti ratkaistua (Subscript[x, i+1]). Kitkatermin kanssa on käytännössä turvauduttava numeeriseen yhtälön ratkaisuun, esim. FindRoot-funktioon.
	(((Subscript[x, i+1])-2(Subscript[x, i])+(Subscript[x, i-1]))/(h^2))=G M ((Subscript[x, i])/(((((Subscript[x, i])^2)+((Subscript[y, i])^2)))^(3/2)))+(γ/m)(Sqrt[((((((Subscript[x, i+1])-(Subscript[x, i-1]))/(2h))))^2)+((((((Subscript[y, i+1])-(Subscript[y, i-1]))/(2h))))^2)])(( (Subscript[x, i+1])-(Subscript[x, i-1]))/(2h))
	(((Subscript[y, i+1])-2(Subscript[y, i])+(Subscript[y, i-1]))/(h^2))=G M ((Subscript[y, i])/(((((Subscript[x, i])^2)+((Subscript[y, i])^2)))^(3/2)))+(γ/m)(Sqrt[((((((Subscript[x, i+1])-(Subscript[x, i-1]))/(2h))))^2)+((((((Subscript[y, i+1])-(Subscript[y, i-1]))/(2h))))^2)])(( (Subscript[y, i+1])-(Subscript[y, i-1]))/(2h))

Ratkaistaan nämä ensin ilman kitkatermejä. Alkunopeus ja -paikka:

θ=60.0Degree;
r0=R{Cos[θ],Sin[θ]};
V=5000;
α=20.0Degree;
v0=V{Cos[α],Sin[α]};
g0=Show[Graphics[Circle[{0,0},R]],Graphics[Point[r0]],Graphics[Arrow[{r0,r0+200v0}]],Axes->True];

Ajanhetket, joilla ratkaisua etsitään:

h=2;t[1]=0.0;T=2000;
i=1;
While[t[i]≤T,i++;t[i]=t[i-1]+h];
n=i;

Alkuehdot ratkaistaville funktioille:

Clear[x,y]
x[1]=r0[[1]];y[1]=r0[[2]];
x[2]=x[1]+v0[[1]]h;
y[2]=y[1]+v0[[2]]h;

Ratkaisuprosessi voidaan keskeyttää, jos kappale päätyy Maan sisään.

Do[
{x[i+1],y[i+1]}=
{x[i+1],y[i+1]}/.
FindRoot[{((x[i+1]-2x[i]+x[i-1])/(h^2))+((G M x[i])/(((((x[i])^2)+((y[i])^2)))^(3/2)))==0,
((y[i+1]-2y[i]+y[i-1])/(h^2))+((G M y[i])/(((((x[i])^2)+((y[i])^2)))^(3/2)))==0},
{x[i+1],x[i]},{y[i+1],y[i]}];
If[Norm[{x[i],y[i]}]<R,Break[]],
{i,2,n-1}]//Quiet;

Show[g0,ListPlot[Table[{x[i],y[i]},{i,1,n}]]]