Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Chapter 06 - Olvaldo corrections #198

Merged
merged 4 commits into from
Apr 10, 2023
Merged

Chapter 06 - Olvaldo corrections #198

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
35a5f76
Apply some Osvaldo corrections
entropidelic Mar 9, 2023
7861061
Apply more corrections
entropidelic Mar 9, 2023
162e4f1
Apply some more corrections
entropidelic Mar 9, 2023
97237e7
Refactor some text
entropidelic Mar 23, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
99 changes: 52 additions & 47 deletions 06_gravity.Rmd
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
# Escaping from Mars

```{julia chap_6_libraries, cache = TRUE, results = FALSE, echo=FALSE}
```{julia chap_6_libraries, cache = FALSE, results = FALSE, echo=FALSE}
cd("./06_gravity/")
import Pkg
Pkg.activate(".")
Expand All @@ -27,32 +27,33 @@ plot(im_0,axis=nothing,border=:none)
#img_small = imresize(im_0, ratio=1/3)
```

To simplify, we approximate the escape velocity as:
Assuming that the planet is spherical and considering only the gravitational effect of Mars, the escape velocity can be described as

$v_{escape}=\sqrt{2*g_{planet}*r_{planet}}$

where $r$ is the radius of the planet and $g$ the constant of gravity at the surface. Suppose that we remember from school that the escape velocity
from Earth is $11\frac{km}{s}$ and that the radius of Mars if half of the earth's.
where $r$ is the radius of the planet and $g$ the gravity constant at the surface of Mars. Suppose that we remember from school that the
escape velocity from Earth is $11\frac{km}{s}$ and that the radius of Mars if half of the earth's.

We remember that the gravity of Earth at its surface is $9.8\frac{m}{s^2}$, so all we need to estimate the escape velocity of Mars is the gravity
of the planet at its surface. So we decide to make an experiment and gather some data. But what exactly do you need to measure? Let's see.

We are going to calculate the constant $g_{mars}$ just throwing stones. We are going to explain a bit the equations regarding the experiment.
The topic we need to revisit is Proyectile Motion.

## Proyectile Motion
Gravity pulls us down to Earth, or in our case, to Mars. This means that we have an acceletation, since there is a force. Recalling the newton equation:
## Proyectile Motion
Gravity pulls us down to Earth, or in our case, to Mars. This means that we have an acceletation, since there is a force. Recalling the
newton equation:

$\overrightarrow{F} = m * \overrightarrow{a}$

where $m$ is the mass of the object, $\overrightarrow{F}$ is the force(it's what make us fall) and $\overrightarrow{a}$ is the acceleration, in
our case is what we call gravity $\overrightarrow{g}$.
The arrow $\overrightarrow{}$ over the letter means that the quantity has a direction in space, in our case, gravity is pointing to the center of
the Earth, or Mars.
where $m$ is the mass of the object, $\overrightarrow{F}$ is the force(it's what make us fall) and $\overrightarrow{a}$ is the acceleration, in
our case is what we call gravity $\overrightarrow{g}$.
The arrow $\overrightarrow{}$ over the letter means that the quantity has a direction in space, in our case, gravity is pointing to the center of
the Earth, or Mars.

*How can we derive the motion of the stones with that equation?*

In the figure below we show a sketch of the problem: We have the 2 axis, $x$ and $y$, the $x$ normally is parallel to the ground and the $y$ axis
In the figure below we show a sketch of the problem: We have the 2 axis, $x$ and $y$, the $x$ axis is parallel to the ground and the $y$ axis
is perpendicular, pointing to the sky. We also draw the initial velocity $v_0$ of the proyectile, and the angle $\theta$ with respect to the ground.
Also it's important to notice that the gravity points in the opposite direction of the $y$ axis.

Expand Down Expand Up @@ -81,14 +82,16 @@ can go. Then, the stone starts to go down.
*How does the velocity evolve in this trajectory?*

Since the begining, the velocity starts decreasing until it has the value of 0 at the highest point, where the stone stops for a moment, then it
changes its direction and start to increase again, pointing towards the ground. Ploting the evolution of the height of the stone, we obtain the plot shown below. We see that, at the begining the stone starts to go up fast and then it slows down. We see that for each value of $y$ there are 2 values of $t$ that satisfies the equation, thats because the stone pass twice for each point, except for the highest value of $y$.
changes its direction and start to increase again, pointing towards the ground. Ploting the evolution of the height of the stone, we obtain
the plot shown below. We see that, at the begining the stone starts to go up fast and then it slows down. We see that for each value of $y$
there are 2 values of $t$ that satisfies the equation, thats because the stone pass twice for each point, except for the highest value of $y$.

```{julia chap_6_plot_3, echo=FALSE}
im_3 = load("./06_gravity/images/img90_.png");
plot(im_3,axis=nothing,border=:none)
```

So, in the example we ahve just explained, we have that the throwing angle is θ=90°, so sin(90°)=1, the trajectory in $y$ becomes:
So, in the example we have just explained, we have that the throwing angle is θ=90°, so sin(90°)=1, the trajectory in $y$ becomes:

$y(t) = v_0*t -\frac{g*t^2}{2}$
And the velocity, which is the derivative of the above equation becomes:
Expand All @@ -108,7 +111,7 @@ The experiment set up will go like this:
- One person will be throwing stones with an angle.
- The other person will be far, watching from some distance, measuring the time since the other throw the stone and it hits the ground. The other
measurement we will need is the distance Δx the stone travelled.
- Also, for the first iteration of the experiment, suppose we only keep the measurements with and initial angle θ~45° (we will loosen this
- Also, for the first iteration of the experiment, suppose we only keep the measurements with an initial angle θ~45° (we will loosen this
constrain in a bit).

```{julia chap_6_plot_4, echo=FALSE}
Expand All @@ -127,14 +130,13 @@ using Turing
Suppose we did the experiment and we have measured then the 5 points, Δx and Δt, shown below:

```{julia, results = FALSE}
Δx_measured = [25.94, 38.84, 52.81, 45.54, 17.24]
t_measured = [3.91, 4.57, 5.43, 4.85, 3.15]
Δx = [25.94, 38.84, 52.81, 45.54, 17.24]
Δt = [3.91, 4.57, 5.43, 4.85, 3.15]
```

*Now, how do we estimate the constant g from those points?*

Using the equations of the trajectory, when the stone hits the ground, $y(t) = 0$, since we take the start of the $y$ coordinate in the ground
(negleting the initial height with respect to the maximum height), so finding the other then the initial point that fulfill this equation, we find that:
When the stone hits the ground, we have that $y(t) = 0$. If we solve for $t$, ignoring the trivial solution $t = 0$, we find that

$t_{f} = \frac{2*v_{0}*sin(θ)}{g}$

Expand All @@ -145,7 +147,7 @@ $Δx=t_{f}*v_{0}*cos(θ)$

where Δx is the distance traveled by the stone.

So, solving for $v_{0}$m, the initial velocity, an unknown quantity, we have:
So, solving for $v_{0}$, the initial velocity, an unknown quantity, we have:

$v_{0}=\frac{Δx}{t_{f}cos(θ)}$

Expand Down Expand Up @@ -188,7 +190,7 @@ values of 0 and 10
```{julia chap_6_plot_6}
plot(Uniform(0,10),xlim=(-1,11), ylim=(0,0.2), legend=false, fill=(0, .5, :lightblue));
xlabel!("g_mars");
ylabel!("Probability");
ylabel!("Probability density");
title!("Uniform prior distribution for g")
```

Expand All @@ -209,21 +211,17 @@ end

```{julia,results = FALSE}
iterations = 10000
ϵ = 0.05
τ = 10
```

```{julia,results = FALSE}
θ = 45
chain_uniform = sample(gravity_uniform(t_measured, Δx_measured, θ), HMC(ϵ, τ), iterations, progress=false);
chain_uniform = sample(gravity_uniform(Δt, Δx, θ), NUTS(), iterations, progress=false);
```

Plotting the posterior distribution for p we that the values are mostly between 2 and 5, with the maximun near 3,8. Can we narrow the values we obtain?
Plotting the posterior distribution for p we that the values are mostly between 2 and 5, with the maximun near 3,8. Can we narrow the values
we obtain?

```{julia chap_6_plot_7}
histogram(chain_uniform[:g], xlim=[1,6], legend=false, normalized=true);
xlabel!("g_mars");
ylabel!("Probability");
ylabel!("Probability density");
title!("Posterior distribution for g with Uniform prior")
```

Expand All @@ -233,7 +231,7 @@ variance of 2, and let the model update its beliefs with the points we have.
```{julia chap_6_plot_8}
plot(Normal(5,2), legend=false, fill=(0, .5,:lightblue));
xlabel!("g_mars");
ylabel!("Probability");
ylabel!("Probability density");
title!("Normal prior distribution for g")
```
We define then the model with a gaussian distribution as a prior for $g$:
Expand All @@ -254,26 +252,27 @@ end
Now we sample values from the posterior distribution and plot and histogram with the values obtained:

```{julia,results = FALSE}
chain_normal = sample(gravity_normal(t_measured, Δx_measured, θ), HMC(ϵ, τ), iterations, progress=false)
chain_normal = sample(gravity_normal(Δt, Δx, θ), NUTS(), iterations, progress=false)
```

```{julia chap_6_plot_9}
histogram(chain_normal[:g], xlim=[3,4.5], legend=false, normalized=true, title="Posterior distribution for g with Normal prior");
xlabel!("g_mars");
ylabel!("Probability");
ylabel!("Probability density");
title!("Posterior distribution for g with Normal prior")
```

We see that the plausible values for the gravity have a clear center in 3.7 and now the distribution is narrower, that's good, but we can do better.
If we observe the prior distribution proposed for $g$, we see that some values are negative, which makes no sense because if that would the
case when you trow the stone, it would go up and up, escaping from the planet.

If we observe the prior distribution proposed for $g$, we see that some values are negative, which makes no sense because if that would the case when you trow the stone, it would go up and up, escaping from the planet.

We propose then a new model for not allowing the negative values to happen. The distribution we are interested in is a LogNormal distribution. In the plot below is the prior distribution for g, a LogNormal distribution with mean 1.5 and variance of 0.5.
We propose then a new model for not allowing the negative values to happen. The distribution we are interested in is a LogNormal distribution.
In the plot below is the prior distribution for g, a LogNormal distribution with mean 1.5 and variance of 0.5.

```{julia chap_6_plot_10}
plot(LogNormal(1,0.5), xlim=(0,10), legend=false, fill=(0, .5,:lightblue));
xlabel!("g_mars");
ylabel!("Probability");
ylabel!("Probability density");
title!("Prior LogNormal distribution for g")
```

Expand All @@ -293,26 +292,28 @@ end
```

```{julia,results = FALSE}
chain_lognormal = sample(gravity_lognormal(t_measured, Δx_measured, θ), HMC(ϵ, τ), iterations, progress=false)
chain_lognormal = sample(gravity_lognormal(Δt, Δx, θ), NUTS(), iterations, progress=false)
```

```{julia chap_6_plot_11}
histogram(chain_lognormal[:g], xlim=[3,4.5], legend=false, title="Posterior distribution for g with LogNormal prior", normalized=true);
xlabel!("g_mars");
ylabel!("Probability");
ylabel!("Probability density");
title!("Posterior distribution for g with LogNormal prior")
```

## Optimizing the throwing angle
Now that we have a good understanding of the equations and the overall problem, we are going to add some difficulties and we will loosen a constrain we have imposed: Suppose that the device employed to measure the angle has an error of 15°, no matter the angle.
Let us remove the angle constraint and add complexity to the problem. We will consider that the device measuring the angle has a
15 degree error, no matter the angle.

*We want to know what are the most convenient angle to do the experiment and to measure or if it doesn't matter.*
Does a most convenient angle to do the experiment exist? Or it doesn’t matter?
Is there an angle that is most convenient for making the experiment? Or doesn't it even matter?

To do the analysis we need to see how the angle influence the computation of $g$, so solving the equation for $g$ we have:

$g = \frac{2*tg(\theta)*\Delta x}{t^{2}_f}$

We can plot then the tangent of θ, with and error of 15° and see what is its maximum and minimun value:
We can plot then the tangent of θ, with a 15 degree error and see what their maximum and minimum values are:

```{julia, results = FALSE}
angles = 0:0.1:70
Expand Down Expand Up @@ -343,8 +344,10 @@ title!("Percentual error");
vline!([angles[findfirst(x->x==minimum(perc_error), perc_error)]], lw=3, label="Minimum error")
```

So, now we see that the lowest percentual error is obtained when we work in angles near 45°, so we are good to go and we can use the data we measured adding the error in the angle.
We now define the new model, where we include an uncertainty in the angle. We propose an uniform prior for the angle centered at 45°, the angle we think the measurement was done.
So, now we see that the lowest percentual error is obtained when we work in angles near 45°, so we are good to go and we can use the data we
measured adding the error in the angle.
We now define the new model, where we include an uncertainty in the angle. We propose an uniform prior for the angle centered at 45°,
the angle we think the measurement was done.

```{julia, results = FALSE}
@model gravity_angle_uniform(t_final, x_final, θ) = begin
Expand All @@ -362,13 +365,13 @@ end
```

```{julia, results = FALSE}
chain_uniform_angle = sample(gravity_angle_uniform(t_measured, Δx_measured, θ), HMC(ϵ, τ), iterations, progress=false)
chain_uniform_angle = sample(gravity_angle_uniform(Δt, Δx, θ), NUTS(), iterations, progress=false)
```

```{julia chap_6_plot_14}
histogram(chain_uniform_angle[:g], legend=false, normalized=true);
xlabel!("g_mars");
ylabel!("Probability");
ylabel!("Probability density");
title!("Posterior distribution for g, including uncertainty in the angle")
```

Expand Down Expand Up @@ -403,18 +406,20 @@ v = 11 .* sqrt.(chain_uniform_angle[:g] ./ (9.8*2))
histogram(v, legend=false, normalized=true);
title!("Escape velocity from Mars");
xlabel!("Escape Velocity of Mars [km/s]");
ylabel!("Probability")
ylabel!("Probability density")
```

Finally, we obtained the escape velocity scape from Mars.

## Summary
In this chapter we had to find the escape velocity from Mars.
To solve this problem, we first needed to find the gravity of Mars, so we started with a physical description of the problem and concluded that by measuring the distance and time of a rock throw plus some Bayesian analysis we could infer the gravity of Mars.
To solve this problem, we first needed to find the gravity of Mars, so we started with a physical description of the problem and concluded
that by measuring the distance and time of a rock throw plus some Bayesian analysis we could infer the gravity of Mars.

Then we created a simple probabilistic model, with the prior probability set to a uniform distribution and the likelihood to a normal distribution.
We sampled the model and obtained our first posterior probability.
We repeated this process two more times, changing the prior distribution of the model for more accurate ones, first with a normal distribution and then with a logarithmic one.
We repeated this process two more times, changing the prior distribution of the model for more accurate ones, first with a Normal distribution and
then with a LogNormal one.

Finally, we used the gravity we inferred to calculate the escape velocity from Mars.

Expand Down
Loading