Discussion of conditions for conservative forces and summing up our discussion on conservative forces. Discussion of potential surfaces and their interpretations.
Reading suggestion: Taylor sections 4.6, 4.9, 4.10 and 5.1 and 5.2 on harmonic oscillations.
The Earth-Sun problem and energy-conserving algorithms and how to encode in more efficient ways various algorithms for solving the equations of motion (Euler, Euler-Cromer and Velocity Verlet).
Reading suggestions: Taylor section 4.8 and these notes
Working on the Earth-Sun problem and hw 5. Hints on various exercises.
Reading suggestions: Taylor chapters 3 and 4 and these notes.
The concept of line integrals plays an important role in our discussion of energy conservation, our definition of potentials and conservative forces.
Let us remind ourselves of some the basic elements (most of you may have seen this in a calculus course under the general topic of vector fields).
We define a path integration \( C \), that is we integrate from a point \( \boldsymbol{r}_1 \) to a point \( \boldsymbol{r}_2 \). Let us assume that the path \( C \) is represented by an arc length \( s \). In three dimension we have the following representation of \( C \)
$$
\boldsymbol{r}(s)=x(s)\boldsymbol{e}_1+y(s)\boldsymbol{e}_2+z(s)\boldsymbol{e}_3,
$$
then our integral of a function \( f(x,y,z) \) along the path \( C \) is defined as
$$
\int_Cf(x,y,z)ds=\int_a^bf\left(x(s),y(s),z(s)\right)ds,
$$
where the initial and final points are \( a \) and \( b \), respectively.
With the definition of a line integral, we can in turn set up the theorem of independence of integration path.
Let us define \( f(x,y,z) \), \( g(x,y,z) \) and \( h(x,y,z) \) to be functions which are defined and continuous in a domain \( D \) in space. Then a line integral like the above is said to be independent of path in \( D \), if for every pair of endpoints \( a \) and \( b \) in \( D \) the value of the integral is the same for all paths \( C \) in \( D \) starting from a point \( a \) and ending in a point \( b \). The integral depends thus only on the integration limits and not on the path.
An expression of the form
$$
fdx+gdy+hdz,
$$
where \( f \), \( g \) and \( h \) are functions defined in \( D \), is a called a first-order differential form in three variables. The form is said to be exact if it is the differential
$$
du= \frac{\partial u}{\partial x}dx+\frac{\partial u}{\partial y}dy+\frac{\partial u}{\partial z}dz,
$$
of a differentiable function \( u(x,y,z) \) everywhere in \( D \), that is
$$
du=fdx+gdy+hdz.
$$
It is said to be exact if and only if we can then set
$$
f=\frac{\partial u}{\partial x},
$$
and
$$
g=\frac{\partial u}{\partial y},
$$
and
$$
h=\frac{\partial u}{\partial z},
$$
everywhere in the domain \( D \).
In vector language the above means that the differential form
$$
fdx+gdy+hdz,
$$
is exact in \( D \) if and only if the vector function (it could be a force, or velocity, acceleration or other vectors we encounter in this course)
$$
\boldsymbol{F}=f\boldsymbol{e}_1+g\boldsymbol{e}_2+h\boldsymbol{e}_3,
$$
is the gradient of a function \( u(x,y,z) \)
$$
\boldsymbol{v}=\boldsymbol{\nabla}u=\frac{\partial u}{\partial x}\boldsymbol{e}_1+\frac{\partial u}{\partial y}\boldsymbol{e}_2+\frac{\partial u}{\partial z}\boldsymbol{e}_3.
$$
If this is the case, we can state the path independence theorem which states that with functions \( f(x,y,z) \), \( g(x,y,z) \) and \( h(x,y,z) \) that fulfill the above exactness conditions, the line integral
$$
\int_C\left(fdx+gdy+hdz\right),
$$
is independent of path in \( D \) if and only if the differential form under the integral sign is exact in \( D \).
This is the path independence theorem.
We will not give a proof of the theorem. You can find this in any vector analysis chapter in a mathematics textbook.
We note however that the path integral from a point \( p \) to a final point \( q \) is given by
$$
\int_p^q\left(fdx+gdy+hdz\right)=\int_p^q\left(\frac{\partial u}{\partial x}dx+\frac{\partial u}{\partial y}dy+\frac{\partial u}{\partial z}dz\right)=\int_p^qdu.
$$
Assume now that we have a dependence on a variable \( s \) for \( x \), \( y \) and \( z \). We have then
$$
\int_p^qdu=\int_{s_1}^{s_2}\frac{du}{ds}ds = u(x(s),y(s),z(s))\vert_{s=s_1}^{s=s_2}=u(q)-u(p).
$$
This last equation
$$
\int_p^q\left(fdx+gdy+hdz\right)=u(q)-u(p),
$$
is the analogue of the usual formula
$$
\int_a^bf(x)dx=F(x)\vert_a^b=F(b)-F(a),
$$
with \( F'(x)=f(x) \).
We remember that a the work done by a force \( \boldsymbol{F}=f\boldsymbol{e}_1+g\boldsymbol{e}_2+h\boldsymbol{e}_3 \) on a displacemnt \( d\boldsymbol{r} \) is
$$
W=\int_C\boldsymbol{F}d\boldsymbol{r}=\int_C(fdx+gdy+hdz).
$$
From the path independence theorem, we know that this has to result in the difference between the two endpoints only. This is exact if and only if the force is the force \( \boldsymbol{F} \) is the gradient of a scalar function \( u \). We call this scalar function, which depends only the positions \( x,y,z \) for the potential energy \( V(x,y,z)=V(\boldsymbol{r}) \).
We have thus
$$
\boldsymbol{F}(\boldsymbol{r})\propto \boldsymbol{\nabla}V(\boldsymbol{r}),
$$
and we define this as
$$
\boldsymbol{F}(\boldsymbol{r})= -\boldsymbol{\nabla}V(\boldsymbol{r}).
$$
Such a force is called a conservative force. The above expression can be used to demonstrate energy conservation.
Finally we can define the criterion for exactness and independence of path. This theorem states that if \( f(x,y,z) \), \( g(x,y,z) \) and \( h(x,y,z) \) are continuous functions with continuous first partial derivatives in the domain \( D \), then the line integral
$$
\int_C\left(fdx+gdy+hdz\right),
$$
is independent of path in \( D \) when
$$
\frac{\partial h}{\partial y}=\frac{\partial g}{\partial z},
$$
and
$$
\frac{\partial f}{\partial z}=\frac{\partial h}{\partial x},
$$
and
$$
\frac{\partial g}{\partial x}=\frac{\partial f}{\partial y}.
$$
This leads to the curl of \( \boldsymbol{F} \) being zero
$$
\boldsymbol{\nabla}\times\boldsymbol{F}=\boldsymbol{\nabla}\times\left(-\boldsymbol{\nabla}V(\boldsymbol{r})\right)=0!
$$
A conservative force \( \boldsymbol{F} \) is a defined as the partial derivative of a scalar potential which depends only on the position,
$$
\boldsymbol{F}(\boldsymbol{r})= -\boldsymbol{\nabla}V(\boldsymbol{r}).
$$
This leads to conservation of energy and a path independent line integral as long as the curl of the force is zero, that is
$$
\boldsymbol{\nabla}\times\boldsymbol{F}=\boldsymbol{\nabla}\times\left(-\boldsymbol{\nabla}V(\boldsymbol{r})\right)=0.
$$
This is taken from homework 4, exercises 5.
A particle is under the influence of a force \( F=-kx+kx^3/\alpha^2 \), where \( k \) and \( \alpha \) are constants and \( k \) is positive.
Determine \( V(x) \) and discuss the motion. It can be convenient here to make a sketch/plot of the potential as function of \( x \).
We assume that the potential is zero at say \( x=0 \). Integrating the force from zero to \( x \) gives
$$
V(x) = -\int_0^x F(x')dx'=\frac{kx^2}{2}-\frac{kx^4}{4\alpha^2}.
$$
The following code plots the potential. We have chosen values of \( \alpha=k=1.0 \). Feel free to experiment with other values. We plot \( V(x) \) for a domain of \( x\in [-2,2] \).
import numpy as np
import matplotlib.pyplot as plt
import math
x0= -2.0
xn = 2.1
Deltax = 0.1
alpha = 1.0
k = 1.0
#set up arrays
x = np.arange(x0,xn,Deltax)
n = np.size(x)
V = np.zeros(n)
V = 0.5*k*x*x-0.25*k*(x**4)/(alpha*alpha)
plt.plot(x, V)
plt.xlabel("x")
plt.ylabel("V")
plt.show()
From the plot here (with the chosen parameters)
and leave the potential well for good.
What happens when the energy of the particle is \( E=(1/4)k\alpha^2 \)? Hint: what is the maximum value of the potential energy?
From the figure we see that the potential has a minimum at at \( x=0 \) then rises until \( x=\alpha \) before falling off again. The maximum potential, \( V(x\pm \alpha) = k\alpha^2/4 \). If the energy is higher, the particle cannot be contained in the well. The turning points are thus defined by \( x=\pm \alpha \). And from the previous plot you can easily see that this is the case (\( \alpha=1 \) in the abovementioned Python code).
We will now venture into a study of a system which is energy conserving. The aim is to see if we (since it is not possible to solve the general equations analytically) we can develop stable numerical algorithms whose results we can trust!
We solve the equations of motion numerically. We will also compute quantities like the energy numerically.
We start with a simpler case first, the Earth-Sun system in two dimensions only. The gravitational force \( F_G \) on the earth from the sun is
$$
\boldsymbol{F}_G=-\frac{GM_{\odot}M_E}{r^3}\boldsymbol{r},
$$
where \( G \) is the gravitational constant,
$$
M_E=6\times 10^{24}\mathrm{Kg},
$$
the mass of Earth,
$$
M_{\odot}=2\times 10^{30}\mathrm{Kg},
$$
the mass of the Sun and
$$
r=1.5\times 10^{11}\mathrm{m},
$$
is the distance between Earth and the Sun. The latter defines what we call an astronomical unit AU.
From Newton's second law we have then for the \( x \) direction
$$
\frac{d^2x}{dt^2}=-\frac{F_{x}}{M_E},
$$
and
$$
\frac{d^2y}{dt^2}=-\frac{F_{y}}{M_E},
$$
for the \( y \) direction.
Here we will use that \( x=r\cos{(\theta)} \), \( y=r\sin{(\theta)} \) and
$$
r = \sqrt{x^2+y^2}.
$$
We can rewrite
$$
F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x,
$$
and
$$
F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y,
$$
for the \( y \) direction.
We can rewrite these two equations
$$
F_{x}=-\frac{GM_{\odot}M_E}{r^2}\cos{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}x,
$$
and
$$
F_{y}=-\frac{GM_{\odot}M_E}{r^2}\sin{(\theta)}=-\frac{GM_{\odot}M_E}{r^3}y,
$$
as four first-order coupled differential equations
$$
\frac{dv_x}{dt}=-\frac{GM_{\odot}}{r^3}x,
$$
$$
\frac{dx}{dt}=v_x,
$$
$$
\frac{dv_y}{dt}=-\frac{GM_{\odot}}{r^3}y,
$$
$$
\frac{dy}{dt}=v_y.
$$
The four coupled differential equations
$$
\frac{dv_x}{dt}=-\frac{GM_{\odot}}{r^3}x,
$$
$$
\frac{dx}{dt}=v_x,
$$
$$
\frac{dv_y}{dt}=-\frac{GM_{\odot}}{r^3}y,
$$
$$
\frac{dy}{dt}=v_y,
$$
can be turned into dimensionless equations or we can introduce astronomical units with \( 1 \) AU = \( 1.5\times 10^{11} \).
Using the equations from circular motion (with \( r =1\mathrm{AU} \))
$$
\frac{M_E v^2}{r} = F = \frac{GM_{\odot}M_E}{r^2},
$$
we have
$$
GM_{\odot}=v^2r,
$$
and using that the velocity of Earth (assuming circular motion) is \( v = 2\pi r/\mathrm{yr}=2\pi\mathrm{AU}/\mathrm{yr} \), we have
$$
GM_{\odot}= v^2r = 4\pi^2 \frac{(\mathrm{AU})^3}{\mathrm{yr}^2}.
$$
The four coupled differential equations can then be discretized using Euler's method as (with step length \( h \))
$$
v_{x,i+1}=v_{x,i}-h\frac{4\pi^2}{r_i^3}x_i,
$$
$$
x_{i+1}=x_i+hv_{x,i},
$$
$$
v_{y,i+1}=v_{y,i}-h\frac{4\pi^2}{r_i^3}y_i,
$$
$$
y_{i+1}=y_i+hv_{y,i},
$$
The code here implements Euler's method for the Earth-Sun system using a more compact way of representing the vectors. Alternatively, you could have spelled out all the variables \( v_x \), \( v_y \), \( x \) and \( y \) as one-dimensional arrays.
# Common imports
import numpy as np
import pandas as pd
from math import *
import matplotlib.pyplot as plt
import os
# Where to save the figures and data files
PROJECT_ROOT_DIR = "Results"
FIGURE_ID = "Results/FigureFiles"
DATA_ID = "DataFiles/"
if not os.path.exists(PROJECT_ROOT_DIR):
os.mkdir(PROJECT_ROOT_DIR)
if not os.path.exists(FIGURE_ID):
os.makedirs(FIGURE_ID)
if not os.path.exists(DATA_ID):
os.makedirs(DATA_ID)
def image_path(fig_id):
return os.path.join(FIGURE_ID, fig_id)
def data_path(dat_id):
return os.path.join(DATA_ID, dat_id)
def save_fig(fig_id):
plt.savefig(image_path(fig_id) + ".png", format='png')
DeltaT = 0.001
#set up arrays
tfinal = 10 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, a, v, and x
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
# Initial conditions as compact 2-dimensional arrays
r0 = np.array([1.0,0.0])
v0 = np.array([0.0,2*pi])
r[0] = r0
v[0] = v0
Fourpi2 = 4*pi*pi
# Start integrating using Euler's method
for i in range(n-1):
# Set up the acceleration
# Here you could have defined your own function for this
rabs = sqrt(sum(r[i]*r[i]))
a = -Fourpi2*r[i]/(rabs**3)
# update velocity, time and position using Euler's forward method
v[i+1] = v[i] + DeltaT*a
r[i+1] = r[i] + DeltaT*v[i]
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots()
#ax.set_xlim(0, tfinal)
ax.set_ylabel('y[AU]')
ax.set_xlabel('x[AU]')
ax.plot(r[:,0], r[:,1])
fig.tight_layout()
save_fig("EarthSunEuler")
plt.show()
We notice here that Euler's method doesn't give a stable orbit. It means that we cannot trust Euler's method. In a deeper way, as we will see in homework 5, Euler's method does not conserve energy. It is an example of an integrator which is not symplectic.
Here we present thus two methods, which with simple changes allow us to avoid these pitfalls. The simplest possible extension is the so-called Euler-Cromer method. The changes we need to make to our code are indeed marginal here. We need simply to replace
r[i+1] = r[i] + DeltaT*v[i]
in the above code with the velocity at the new time \( t_{i+1} \)
r[i+1] = r[i] + DeltaT*v[i+1]
By this simple caveat we get stable orbits. Below we derive the Euler-Cromer method as well as one of the most utlized algorithms for sovling the above type of problems, the so-called Velocity-Verlet method.
Let us repeat Euler's method. We have a differential equation
$$
\begin{equation}
y'(t_i)=f(t_i,y_i)
\tag{1}
\end{equation}
$$
and if we truncate at the first derivative, we have from the Taylor expansion
$$
\begin{equation}
y_{i+1}=y(t_i) + (\Delta t) f(t_i,y_i) + O(\Delta t^2), \tag{2}
\end{equation}
$$
which when complemented with \( t_{i+1}=t_i+\Delta t \) forms the algorithm for the well-known Euler method. Note that at every step we make an approximation error of the order of \( O(\Delta t^2) \), however the total error is the sum over all steps \( N=(b-a)/(\Delta t) \) for \( t\in [a,b] \), yielding thus a global error which goes like \( NO(\Delta t^2)\approx O(\Delta t) \).
To make Euler's method more precise we can obviously decrease \( \Delta t \) (increase \( N \)), but this can lead to loss of numerical precision. Euler's method is not recommended for precision calculation, although it is handy to use in order to get a first view on how a solution may look like.
Euler's method is asymmetric in time, since it uses information about the derivative at the beginning of the time interval. This means that we evaluate the position at \( y_1 \) using the velocity at \( v_0 \). A simple variation is to determine \( x_{n+1} \) using the velocity at \( v_{n+1} \), that is (in a slightly more generalized form)
$$
\begin{equation}
y_{n+1}=y_{n}+ v_{n+1}+O(\Delta t^2)
\tag{3}
\end{equation}
$$
and
$$
\begin{equation}
v_{n+1}=v_{n}+(\Delta t) a_{n}+O(\Delta t^2).
\tag{4}
\end{equation}
$$
The acceleration \( a_n \) is a function of \( a_n(y_n, v_n, t_n) \) and needs to be evaluated as well. This is the Euler-Cromer method.
Exercise: go back to the above code with Euler's method and add the Euler-Cromer method.
The Euler-Cromer method conserves energy. If you are interested, the following article shows why.
Let us stay with \( x \) (position) and \( v \) (velocity) as the quantities we are interested in.
We have the Taylor expansion for the position given by
$$
x_{i+1} = x_i+(\Delta t)v_i+\frac{(\Delta t)^2}{2}a_i+O((\Delta t)^3).
$$
The corresponding expansion for the velocity is
$$
v_{i+1} = v_i+(\Delta t)a_i+\frac{(\Delta t)^2}{2}v^{(2)}_i+O((\Delta t)^3).
$$
Via Newton's second law we have normally an analytical expression for the derivative of the velocity, namely
$$
a_i= \frac{d^2 x}{dt^2}\vert_{i}=\frac{d v}{dt}\vert_{i}= \frac{F(x_i,v_i,t_i)}{m}.
$$
If we add to this the corresponding expansion for the derivative of the velocity
$$
v^{(1)}_{i+1} = a_{i+1}= a_i+(\Delta t)v^{(2)}_i+O((\Delta t)^2)=a_i+(\Delta t)v^{(2)}_i+O((\Delta t)^2),
$$
and retain only terms up to the second derivative of the velocity since our error goes as \( O(h^3) \), we have
$$
(\Delta t)v^{(2)}_i\approx a_{i+1}-a_i.
$$
We can then rewrite the Taylor expansion for the velocity as
$$
v_{i+1} = v_i+\frac{(\Delta t)}{2}\left( a_{i+1}+a_{i}\right)+O((\Delta t)^3).
$$
Our final equations for the position and the velocity become then
$$
x_{i+1} = x_i+(\Delta t)v_i+\frac{(\Delta t)^2}{2}a_{i}+O((\Delta t)^3),
$$
and
$$
v_{i+1} = v_i+\frac{(\Delta t)}{2}\left(a_{i+1}+a_{i}\right)+O((\Delta t)^3).
$$
Note well that the term \( a_{i+1} \) depends on the position at \( x_{i+1} \). This means that you need to calculate the position at the updated time \( t_{i+1} \) before the computing the next velocity. Note also that the derivative of the velocity at the time \( t_i \) used in the updating of the position can be reused in the calculation of the velocity update as well.
We can now easily add the Verlet method to our original code as
DeltaT = 0.01
#set up arrays
tfinal = 10 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, a, v, and x
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
# Initial conditions as compact 2-dimensional arrays
r0 = np.array([1.0,0.0])
v0 = np.array([0.0,2*pi])
r[0] = r0
v[0] = v0
Fourpi2 = 4*pi*pi
# Start integrating using the Velocity-Verlet method
for i in range(n-1):
# Set up forces, air resistance FD, note now that we need the norm of the vecto
# Here you could have defined your own function for this
rabs = sqrt(sum(r[i]*r[i]))
a = -Fourpi2*r[i]/(rabs**3)
# update velocity, time and position using the Velocity-Verlet method
r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
rabs = sqrt(sum(r[i+1]*r[i+1]))
anew = -4*(pi**2)*r[i+1]/(rabs**3)
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
t[i+1] = t[i] + DeltaT
# Plot position as function of time
fig, ax = plt.subplots()
ax.set_ylabel('y[AU]')
ax.set_xlabel('x[AU]')
ax.plot(r[:,0], r[:,1])
fig.tight_layout()
save_fig("EarthSunVV")
plt.show()
You can easily generalize the calculation of the forces by defining a function which takes in as input the various variables. We leave this as a challenge to you.
Here we show how we can set up a code which conserves the total energy. We need then to compute the kinetic energy and the potential energy. Here we have defined two simple arrays as functions of time for the kinetic and the potential energies which we update as soon as we have updated the positions and the velocities. Note that we are keeping the Sun at rest, meaning that it is only the Earth which moves.
DeltaT = 0.0001
#set up arrays
tfinal = 10 # in years
n = ceil(tfinal/DeltaT)
# set up arrays for t, a, v, and x
t = np.zeros(n)
v = np.zeros((n,2))
r = np.zeros((n,2))
Ekin = np.zeros(n)
Epot = np.zeros(n)
# Initial conditions as compact 2-dimensional arrays
r0 = np.array([1.0,0.0])
v0 = np.array([0.0,2.0*pi])
r[0] = r0
v[0] = v0
Fourpi2 = 4*pi*pi
# Mass of Earth set to 1 and setting up initial Ekin and Epot
Ekin[0] = 0.5*sum(v[0]*v[0])
Epot[0] = -Fourpi2/sqrt(sum(r[0]*r[0]))
# Start integrating using the Velocity-Verlet method
for i in range(n-1):
# Set up forces, air resistance FD, note now that we need the norm of the vecto
# Here you could have defined your own function for this
rabs = sqrt(sum(r[i]*r[i]))
a = -Fourpi2*r[i]/(rabs**3)
# update velocity, time and position using the Velocity-Verlet method
r[i+1] = r[i] + DeltaT*v[i]+0.5*(DeltaT**2)*a
rabs = sqrt(sum(r[i+1]*r[i+1]))
anew = -4*(pi**2)*r[i+1]/(rabs**3)
v[i+1] = v[i] + 0.5*DeltaT*(a+anew)
Ekin[i+1] += 0.5*(sum(v[i+1]*v[i+1]))
Epot[i+1] += -Fourpi2/rabs
t[i+1] = t[i] + DeltaT
# Plot energy as function of time
fig, ax = plt.subplots()
ax.set_ylabel('E')
ax.set_xlabel('t[yr]')
ax.plot(t, Epot+Ekin)
fig.tight_layout()
save_fig("EnergyEarthSunVV")
plt.show()
This exercise is discussed in Taylor's chapter 3.2.
Consider the rocket of mass \( M \) moving with velocity \( v \). After a brief instant, the velocity of the rocket is \( v+\Delta v \) and the mass is \( M-\Delta M \). Momentum conservation gives
$$
\begin{eqnarray*}
Mv&=&(M-\Delta M)(v+\Delta v)+\Delta M(v-v_e)\\
0&=&-\Delta Mv+M\Delta v+\Delta M(v-v_e),\\
0&=&M\Delta v-\Delta Mv_e.
\end{eqnarray*}
$$
In the second step we ignored the term \( \Delta M\Delta v \) since we can assume it is small. The last equation gives
$$
\begin{eqnarray}
\Delta v&=&\frac{v_e}{M}\Delta M,\\
\nonumber
\frac{dv}{dt}dt&=&\frac{v_e}{M}dM.
\end{eqnarray}
$$
Here we let \( \Delta v\rightarrow dv \) and \( \Delta M\rightarrow dM \). We have also assumed that \( M(t) = M_0-kt \). Integrating the expression with lower limits \( v_0=0 \) and \( M_0 \), one finds
$$
\begin{eqnarray*}
v&=&v_e\int_{M_0}^M \frac{dM'}{M'}\\
v&=&v_e\ln(M/M_0)\\
&=&v_e\ln[(M_0-k t)/M_0].
\end{eqnarray*}
$$
We have ignored gravity here. If we add gravity as the external force, we get when integrating an additional terms \( -gt \), that is
$$
v(t)=v_e\ln[(M_0-k t)/M_0]-gt.
$$
This is a continuation of the previous exercise and most of the relevant background material can be found in Taylor chapter 3.2.
Taking the velocity from the previous exercise and integrating over time we find the height
$$
y(t) = y(t_0=0)+\int_0^tv(t')dt'.
$$
You need to insert \( v(t) \) from the previous exercise
To do the integral over time we recall that \( M(t)=M_0-\Delta M t \). We assumed that \( \Delta M=k \) is a constant. We use that \( M_0-M=kt \) and assume that mass decreases by a constant \( k \) times time \( t \).
We will need to compute an integral which goes like
$$
\int_0^t \ln{M(t')}dt' = \int_0^t \ln{(M_0-kt')}dt'.
$$
and defining the variable \( u=M_0-kt' \), with \( du=-kdt' \) and the new limits \( M_0 \) when \( t=0 \) and \( M_0-kt \) when time is equal to \( t \), we have
$$
\int_0^t \ln{M(t')}dt' = \int_0^t \ln{(M_0-kt')}dt'=-\frac{1}{k}\int_{M_0}^{M_0-kt} \ln{(u)}du=-\frac{1}{k}\left[u\ln{(u)}-u\right]_{M_0}^{M_0-kt}.
$$