# Announcements
- No Problem Set this week, Problem Set 4 will be posted on 9/28.
- Solutions to Problem Set 3 posted on D2L.

<style>
@import url(https://www.numfys.net/static/css/nbstyle.css);
</style>
<a href="https://www.numfys.net"></a>

# Ordinary Differential Equations - higher order methods in practice

<section class="post-meta">
Based on notes and notebooks by Niels Henrik Aase, Thorvald Ballestad, Vasilis Paschalidis and Jon Andreas Støvneng
</section>


### Higher Order Derivatives and Sets of 1st order ODEs

The trick to solving ODEs with higher derivatives is turning them into systems of first-order ODEs.

As a simple example, consider the second-order differential equation describing the van der Pol oscillator

$$ \frac{d^2 x}{dt^2} - a (1-x^2) \frac{dx}{dt} + x = 0 $$

We turn this into a pair of first-order ODEs by defining an auxiliary function $v(t) = dx/dt$ and writing the system as

\begin{align}
\begin{split}
\frac{dv}{dt} &= a (1-x^2) v - x\\
\frac{dx}{dt} &= v
\end{split}
\end{align}

Note that there are only functions (and the independent variable) on the RHS; all "differentials" are on the LHS.

Now that we have a system of first-order equations ,we can proceed as above. A function describing
the RHS of this system is


```python
import numpy as np
import matplotlib.pyplot as plt
def RK4_step(t, y, h, g, *P):
    """
    Implements a single step of a fourth-order, explicit Runge-Kutta scheme
    """
    thalf = t + 0.5*h
    k1 = h * g(t, y, *P)
    k2 = h * g(thalf, y + 0.5*k1, *P)
    k3 = h * g(thalf, y + 0.5*k2, *P)
    k4 = h * g(t + h, y + k3, *P)
    return y + (k1 + 2*k2 + 2*k3 + k4)/6
```


```python
def vdpRHS(t, y, a):
    # we store our function as the array [x, x']
    return np.array([
        
        y[1],                          # dx/dt = v
        a*(1-y[0]**2)*y[1] - y[0]      # dv/dt = a*(1-x**2)*v - x
        
    ])
```


```python
def odeSolve(t0, y0, tmax, h, RHS, method, *P):
    """
    ODE driver with constant step-size, allowing systems of ODE's
    """
    # make array of times and find length of array
    t = np.arange(t0,tmax+h,h)
    ntimes,  = t.shape

    # find out if we are solving a scalar ODE or a system of ODEs, and allocate space accordingly
    
    if type(y0) in [int, float]:  # check if primitive type -- means only one eqn
        neqn = 1
        y = np.zeros( ntimes )
    else:                         # otherwise assume a numpy array -- a system of more than one eqn
        neqn, = y0.shape
        y = np.zeros( (ntimes, neqn) )

    # set first element of solution to initial conditions (possibly a vector) 
    y[0] = y0

    # march on...
    for i in range(0,ntimes-1):
        y[i+1] = method(t[i], y[i], h, RHS, *P)

    return t,y
```


```python
a = 15 # parameter

h = 0.01
t0 = 0.0
y0 = np.array([ 0,  1])
tmax = 3

# call the solver
t, y = odeSolve(t0, y0, tmax, h, vdpRHS, RK4_step, a)

fig,ax = plt.subplots()
ax.plot(t,y[:,0],'b', label='x')
ax.plot(t,y[:,1],'r--', label='v')
ax.set_xlabel('time')
ax.legend()
ax.set_title(f"van der Pol Oscillator for a={a}")

plt.tight_layout()
plt.show()

```

A somewhat more complex example is the Lane-Emden equation, which is really just Poisson's equation in spherical symmetry for the graviational potential of a self-gravitating fluid whose pressure is related to its density as $P\propto\rho^\gamma$. Such a system is called a _polytrope_, and is often used in astrophysics as a simple model for the structure of system such as a a star in which outward pressure and inward gravity are in equilibrium.

Let $\xi$ be the dimensionless radius of the system, and let $\theta$ be related to the density as
$\rho = \rho_c \theta^n$, where $\rho_c$ is the density at the origin and $n = 1/(\gamma-1)$. We then have the dimensionless second-order differential equation

$$ \frac{1}{\xi^2}\frac{d}{d\xi}\left(\xi^2\frac{d\theta}{d\xi}\right) + \theta^n = 0 $$

Note that the first term is just the divergence $\nabla\cdot\theta$ in spherical symmetry.

If we expand out the first term, we have 

$$ \frac{d^2\theta}{d\xi^2} + \frac{2}{\xi}\frac{d\theta}{d\xi} + \theta^n = 0 $$

Defining an auxiliary function $v(\xi) = d\theta/d\xi$, we can then convert this into a system of two first-order ODEs:

\begin{align}
\begin{split}
\frac{dv}{d\xi} &= -\frac{2}{\xi} v - \theta^n \\
\frac{d\theta}{d\xi} &= v
\end{split}
\end{align}

Again, we have "derivatives" only on the LHS and no derivatives on the RHS of our system.

Looking at this expression, one can see right away that at the origin $\xi=0$ we will have a numerical problem; we are dividing by zero.
Analytically, this is not a problem, since $v/\xi\rightarrow0$ as $\xi\rightarrow0$, but here we need to address this numerically.

The first approach is to take care of the problem in our RHS function:


```python
def leRHS(x, y, n):
    
    dthetadx = y[1]
    
    if x==0:
        dvdx = -y[0]**n
    else:
        dvdx = -2/x*y[1] - y[0]**n
    
    return np.array([ dthetadx, dvdx ])

```

This is somewhat clunky, however, and you would first have to convince yourself that in fact $v(\xi)\rightarrow0$ faster than $\xi$ (don't just take my word for it!).

Instead, we could use a more direct RHS function


```python
def leRHS(x, y, n):
    
    dthetadx = y[1]
    
    dvdx = -2/x*y[1] - y[0]**n
    
    return np.array([ dthetadx, dvdx ])
```

and expand the solution in a Taylor series about the origin to get a starting value for our numerical integration at a small distance away from the origin. To do this, write 

$$\theta(\xi) = a_0 + a_1 \xi + a_2 \xi^2 + \dots$$

The first thing to notice is that, by symmetry, only even powers of $\xi$ will appear in the solution.
Thus we will have

$$ \theta(\xi) = a_0 + a_2 \xi^2 + a_4 \xi^4 + \dots$$

By the boundary condition $\theta(0) = 1$, we have immediately that $a_0 = 1$.

Next, substitute $\theta(\xi) = 1 + a_2 \xi^2 + a_4 \xi ^4 + O(\xi^6)$ into the Lane-Emden equation. $\theta$ and its first two derivatives are

\begin{align}
\begin{split}
\theta(\xi) &= 1 + a_2 \xi^2 + a_4 \xi^4 + O(\xi^6)\\
\theta'(\xi) &= 2 a_2 \xi + 4 a_4 \xi^3 + O(\xi^5) \\
\theta''(\xi) &= 2 a_2 + 12 a_4 \xi^2 + O(\xi^4)
\end{split}
\end{align}

Putting these into the Lane-Emden equation, we have
\begin{align}
\begin{split}
2 a_2 + 12 a_4 \xi^2 + O(\xi^4) + \frac{2}{\xi} (2 a_2 x + 4 a_4 \xi^3 + O(\xi^5))  &= -\theta^n \\
6 a_2 + 20 a_4 \xi^2 + O(\xi^4) &= -\theta^n
\end{split}
\end{align}

A boundary condition $\theta(0)=1$, and thus we have $a_2 = -1/6$. Away from zero, then, we have

\begin{align}
\begin{split}
-1 + 20 a_4 \xi^2 + O(\xi^4) &= -\left(1 - 1/6 \xi^2 + a_4 \xi^4 + O(\xi^6)\right)^n
\end{split}
\end{align}

The term on the RHS is $ 1 - n \xi^2/6 + O(\xi^4)$, and so we must have $a_4 = n/120$.

Thus, the series expansion of the solution around the origin is

$$ \theta(\xi) = 1 - \frac{1}{6}\xi^2 + \frac{n}{120} \xi^4 + \dots $$

We can now use this expansion to take a first step slightly away from the origin before beginning our
numerical integration, thus avoiding the divide by zero. Note that this series solution near the origin is $O(h^5)$ and so is a good match for RK4 if we take the same (or smaller) step-size.


```python
n = 3

xi0 = 0.01                            # starting value of xi for our numerical integration
theta0 = 1 - xi0**2/6 + n*xi0**4/120  # Taylor series solution to the DE near zero derived above
theta0p = -xi0/3 + n*xi0**3/30
y0 = np.array([ theta0,  theta0p])    # set IC's for numerical integration
print(f"IC at {xi0:10.5e}: {y0[0]:10.5e}, {y0[1]:10.5e}")

h = 0.1
tmax = 8

# call the solver
t, y = odeSolve(xi0, y0, tmax, h, leRHS, RK4_step, n)

fig,ax = plt.subplots()
ax.plot(t,y[:,0],'b', label=r'$\theta(\xi)$')
ax.plot(t,y[:,1],'r--', label=r'$\frac{d\theta}{d\xi}$')
ax.plot([0,tmax],[0,0],'k')
ax.set_xlabel(r'$\xi$')
ax.set_title(f"Lane Emden Equation for n={n}")
ax.legend()

plt.tight_layout()
plt.show()
```

For values of $n\le5$, the solutions of the Lane Emden equation (the so-called Lane-Emden functions of index $n$) decrease to zero at finite $\xi$. Since this is the radius at which the density goes to zero, we can interpret it as the surface of the self-gravitating body (for example, the radius of the star). Knowing this value $\xi_1$ is thus interesting... Let us see how to determine it numerically.

Cleary, we are looking for the solution to $\theta(\xi_1)=0$; this is just root-finding, which we already know how to do. Instead of using some closed-form function, however, the value of the function $\theta(\xi)$ must in this case be determined numerically. But we have just figured out how to do this!

Let's use the bisection method for our root-finding algorithm; here is a quick version (no error checking!)


```python
def bisection(func, low, high, eps, *P):
    flow = func(low, *P)
    fhigh = func(high, *P)
    mid = 0.5*(low+high)
    fmid = func(mid,*P)
    
    while (high-low)> eps:
        if fmid*flow < 0:
            high = mid
            fhigh = fmid
        else:
            low = mid
            flow = mid
            
        mid = 0.5*(low+high)
        fmid = func(mid,*P)
            
    return low
```

Now let us make a function which returns $\theta(\xi)$, the solution to the Lane-Emden equation at $\xi$


```python
def theta(xi, n):
    h = 1e-4
    xi0 = 1e-4
    theta0 = 1 - xi0**2/6 + n*xi0**4/120
    theta0p = -xi0/3 + n*xi0**3/30
    y0 = np.array([ theta0,  theta0p])
    t, y = odeSolve(xi0, y0, xi, h, leRHS, RK4_step, n)
    
    return y[-1,0]
```

Using these, we can compute the surface radius of the polytrope


```python
n = 3
xi1 = bisection(theta, 6, 8, 1e-5, n)
print(f"xi_1 = {xi1:7.5f}")
```

    xi_1 = 6.89680


A more careful treatment gives a value $\xi_1 = 6.89685...$, so we are doing pretty well...

# Breaktout exercise: Projectile motion

The trajectory of a projectile in a uniform gravity field is a classic and dear physics problem. However, in most introductory physics courses, one works under the assumption that the frictional forces can be neglected. This is of course not always the case, especially not at high speeds. If frictional forces is included, the path of the projectile can no longer be described analytically, and the trajectory must be calculated numerically. To do this, one must solve a coupled set of differential equations, which provides us with an excellent opportunity to implement an ODE-solver on vector form. This turns out to be an elegant and easy implementation, that can be used for solving both differential equations of higher order and coupled differential equations. 

## Theory
### Equations of motion
The governing equations for a projectile is derived from Newtons second law, $\vec{F} = m\vec{a}$, where $\vec{F}$ is the force acting on the projectile, $m$ is the mass, and $\vec{a}$ is the acceleration of the ball. In this example there are two forces acting on the ball. The gravitational force 

$$\vec{F_G} = m\vec{g},$$

where $\vec{g}$ denotes the the standard gravity. The other force is the drag force

$$\vec{F_D} = -\frac{1}{2}\rho ACv_r^2\hat{v_r}.$$ 

$A$ is the cross-sectional area, $C$ is the drag coefficient, and $v_r$ is the velocity of the projectile relative to the fluid it is traveling through, $\rho$ is the air density and $\hat{v_r}$ is the direction of the projectile with length unity. A more practical way of writing this equation is

$$\vec{F_D} = - \frac{B}{\rho_0}\rho v_r\left(\vec{v}-\vec{V}\right),$$

where $\rho_0$ is the air density at sea level altitude, $B$ is some constant, $\vec{v}$ is the velocity of the projectile, and $\vec{V}$ is the velocity of the wind. 

In this example we will use the shells from the German cannon "the Paris gun" as the projectiles, and the numerical constants mentioned above will be based on historical data of the cannon. The Paris gun operated during World War I and bombarded Paris from far away, the distances often exceeded 100 kilometers, with the velocities of the shells being correspondingly high.

The complete equation for the drag force includes an extra term, which is not mentioned above, namely the term originating from Stokes' law. For a spherical projectile the term takes the form 

$$\vec{F_D'} = -6 \pi \eta \vec{v},$$

where $\vec{F_D'}$ denotes the extra term, and $\eta$ denotes the air's viscosity. At high velocities this term is dwarfed by $\vec{F_D}$ as $\vec{F_D}$ is proportional to $v^2$, thus making the drag force $\vec{F_D}$ an excellent approximation for the total drag force. 

The last part of the theory needed before we can find the equations of motion is an expression for the air density. As the shells of the Paris gun will reach the stratosphere, it is not reasonable to treat $\rho$ as a constant. In this notebook we will use an adiabatic air density model. The air density is then given by 

$$\rho = \rho_0 \left(1 - \frac{ay}{T_0}\right)^\alpha,$$

where $\rho$ and $T_0$ denotes the air density and temperature at sea level, $y$ denotes the altitude, $a$ and $\alpha$ are constants $a= 6.5 \cdot 10^{-3}$ K/m, $\alpha$ = 2.5. At high altitudes the expression in the parentheses will become negative, thus making the air density an imaginary number. To avoid this, the air density is set equal to zero when reaching these altitudes. This critical height is around 44 kilometers, well into the stratosphere. 

In a two-dimensional coordinate system, where $x$ and $y$ are the displacement from the origin in horizontal and vertical direction respectively, the equations of motion become
$$
\begin{align}
m\ddot{x} &= -\frac{B}{\rho_0}\rho\left(\dot{x}-V_x\right)\left|\vec{v}-\vec{V}\right|,\\
m\ddot{y} &= -\frac{B}{\rho_0}\rho\left(\dot{y}-V_y\right)\left|\vec{v}-\vec{V}\right|-mg,
\end{align}
$$



where $V_x$ and $V_y$ denotes the velocity of the wind in the horizontal and vertical direction. Simplifying even further, and introducing $u$ and $v$ as the velocities in the horizontal and vertical direction yields four coupled differential equations

\begin{equation}
\begin{aligned}
\dot{x} &= u, \\
\dot{u} &= -\frac{B}{m}\Big[ 1 - \frac{ay}{T_0} \Big]^\alpha \left(u-V_x\right)\Big[ (u-V_x)^2+(v-V_y)^2 \Big]^{1/2},\\
\dot{y} &= w, \\
\dot{w} &= -\frac{B}{m}\Big[ 1 - \frac{ay}{T_0} \Big]^\alpha \left(v-V_y\right)\Big[ (u-V_x)^2+(v-V_y)^2 \Big]^{1/2}-g.\\
\end{aligned}
\label{ODEs}
\end{equation}

### Discretization: Runge-Kutta 4th order

First we define a set of $N$ discrete time values $$
t_n=t_0+n\cdot h~~~~~\mathrm{with}~~~~n=0,1,2,3,...,N.
$$

$h$ is the time step length and is calculated by the formula $h = \frac{t_N-t_0}{N}$, where $t_0$ and $t_N$ denotes the initial and final time values of the projectile's flight. These $N$ time values correspond to $N$ number of $x$, $y$, $u$, and $w$ values, each denoted $x_n$, $y_n$, $u_n$ and $w_n$. A more practical way of denoting this collection of values is on vector form
$$
\vec{w_n} = \begin{bmatrix}x_n \\ y_n \\ u_n \\ v_n\end{bmatrix}.
$$

This notation makes it possible to implement the Runge-Kutta method in an elegant way. It is always possible to reduce an ordinary differential equation of order $q$ to a problem of solving a set of two $q-1$ order coupled differential equations. We have a coupled set of two second order differential equations, thus we need to solve a set of four first order coupled differential equations. By introducing a function $f$ that transforms $\vec{w_n}$ in the way described by equation \eqref{ODEs}, one can implement the Runge-Kutta method directly. In our case, $f$ is given by

$$
f(\vec{w_n} ) = f \begin{bmatrix}x_n \\ y_n \\ u_n \\ v_n\end{bmatrix} = \begin{bmatrix}u_n \\ v_n \\ -\frac{B}{m}\Big[ 1 - \frac{ay_n}{T_0} \Big]^\alpha \left(u_n-V_x\right)\Big[ (u_n-V_x)^2+(v_n-V_y)^2 \Big]^{1/2} \\ -\frac{B}{m}\Big[ 1 - \frac{ay_n}{T_0} \Big]^\alpha \left(v_n-V_y\right)\Big[ (u_n-V_x)^2+(v_n-V_y)^2 \Big]^{1/2}-mg\end{bmatrix},
$$
with no time dependence, such that 
$$
\dot{\vec{w_n}}=f(\vec{w_n})
$$

describes the complete system of differential equations. The following code implements this, and can easily be expanded to solve either an even higher order ODE or an ODE describing a three-dimensional problem. 



```python
# Importing necessary packages
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import rc
%matplotlib inline

newparams = {'figure.figsize': (15, 7), 'axes.grid': False,
             'lines.markersize': 10, 'lines.linewidth': 2,
             'font.size': 15, 'mathtext.fontset': 'stix',
             'font.family': 'STIXGeneral', 'figure.dpi': 200}
plt.rcParams.update(newparams)

# Constants

g = 9.81             # Force of gravity per kilo on earth's surface
alpha = 2.5          # Parameter in the adiabatic air density model
a = 6.5 * 10 ** (-3) # Parameter in the adiabatic air density model
T_0 = 288            # Temperature at sea level
m = 50               # Mass of the projectile
B = 2 * 10 ** (-3)   # Constant based on the Paris gun
v_0 = 1640           # Initial velocity


# Wind strength. This is a constant that never varies 
V = np.zeros(2) # Starting of with no wind


def f(t,w,*P):
    """ A function describing the right hand side of the equations
    of motion/the set of ODEs.
    Parameter:
        t       (time parameter, not needed in this case)
        w       vector containg the needed coordinates and their velocities 
    """
    # V : the vector describing the wind strength and direction
    temp_vector = np.zeros(4)
    temp_vector[0] = w[2]
    temp_vector[1] = w[3]
    # Saving the next parameters in order to optimize run time
    k = (1 - a * w[1] / T_0)   # Air density factor
    s = np.sqrt((w[2]-V[0]) ** 2 + (w[3]-V[1]) ** 2 )
    if k>0:
        temp_vector[2] = - B / m * k ** (alpha) * (w[2]-V[0]) * s
        temp_vector[3] = - B / m * k ** (alpha) * (w[3]-V[1]) * s - g
    else:
        temp_vector[2] = 0 
        temp_vector[3] = - g
    return temp_vector 
```

Now we can finally start to explore the trajectories of the Paris gun. The functions below will be helpful going forward.


```python
def shoot(theta, v0):
    """ Initializes the vector w (x and y position and velocity of the
    projectile) given a initial shooting angle, theta, and 
    absolute velocity, v0. 
    """
    w = np.zeros(4)
    w[2] = v0 * np.cos(np.deg2rad(theta))
    w[3] = v0 * np.sin(np.deg2rad(theta))
    return w


def projectile_motion(h, theta):
    """ Calculates the motion of the projectile using the functions
    defined above. While the projectile is in the air (w[1] >=0) the 
    position and velocity is updated using a single step of RK4.
    Parameters:
        h        unit/integration step length
        theta    initial shooting angle
    Returns:
        X_list   array of the projectile's x-position
        Y_list   array of the porjectile's y-position
    """
    w = shoot(theta, v_0)
    X_list = np.zeros(0)
    Y_list = np.zeros(0)
    while  w[1] >= 0:
        w = RK4_step(0, w, h,f,0)
        X_list = np.append(X_list, w[0])
        Y_list = np.append(Y_list, w[1])
    return X_list, Y_list
```

### Finding the optimal launch angle
With no drag force, it is trivial to show that the angle maximizing the range of a projectile is $\theta = 45^\circ$, where $\theta$ is the angle between the horizontal axis and the direction of the initial velocity. It is interesting to se how the air drag affects this maximum angle, both with and without wind. The following code finds the optimal angle to one decimal precision, and shows how the range depends on the firing angle. It is also instructive to see some of the trajectories of the projectile, so some of the trajectories are included. 


```python
def find_optimal_angle(h):
    """ Given an integration time step, this function calculates the optimal initial 
    shooting angle for the projectile to obtain maximum range, in x-direction. The 
    set of angles tested, with their corresponding range, along with the optimal
    angle are returned.
    """
    
    record = 0        # Placeholder variable that holds the maximum range
    optimal_angle = 0 # Placeholder variable that holds the angle yielding the maximum range
    # Lists containing the initial angle and its corresponding range
    theta_list = np.zeros(0) 
    range_list = np.zeros(0)
    for theta in range (1,90,2):
        x_list, y_list = projectile_motion(h, theta)
        # Using linear interpolation do determine the landing point more precisely
        m = (y_list[-1] - y_list[-2]) / (x_list[-1] - x_list[-2])   # The landing point
        x_range = - y_list[-1] / m + x_list[-1] 
        theta_list = np.append(theta_list, theta)
        range_list = np.append(range_list, x_range)
        # Update records
        if x_range >= record:
            record = x_range
            optimal_angle = theta

    # Rerunning the same code on a smaller interval in order to approximate the optimal angle
    # more precicely    
    return theta_list, range_list, optimal_angle
    
theta, x , best = find_optimal_angle(0.1)

print("The optimal angle is: ", best, " degrees")
```

    The optimal angle is:  55  degrees



```python
plt.plot(theta, x/1000)
plt.title(r"Projectile range as a function of shooting angle, $\theta$")
plt.xlabel(r"$\theta $ [$\degree$]")
plt.ylabel(r"range [km]")

def trajectories(h):
    plt.figure()
    plt.title("Projectile trajectories by alternating shooting angle")
    plt.xlabel(r"x [km]")
    plt.ylabel(r"y [km]")
    theta_list = np.arange(30.0,75,5)
    for angle in theta_list:
        x_list, y_list = projectile_motion(h, angle)
        plt.plot(x_list/1000, y_list/1000, '--', label=r'$\theta = $%.i $\degree$'%(angle))
    plt.legend(loc='best')
    plt.gca().set_ylim(bottom=0)
    plt.show()
trajectories(0.1)
```

In contrast to projectile motion with no friction, maximum range is not obtained with firing angle $\theta = 45^\circ$. Because  the drag force is lower at higher altitudes, the range increases when the projectile spends the majority of it's flight at higher altitudes. This explains why the maximum range is obtained at $\theta = 55^\circ$.

### Numerical validity

In order to determine the validity of the numerical implementation, we can use conservation of energy. During the projectile's trajectory the drag force $\vec{F_D}$ will do work on the projectile. The cumulative work done on the projectile is given by 

$$
W = \int \vec{F_D} \cdot \vec{dr}.
$$

In our discrete model, we can easily calculate this integral by summing over all the infinitesimal contributions, and using the fact that $\vec{F_D}$ always is antiparallel with $\vec{dr}$ such that the cumulative work at time $t_M$, $M<N$, is equal to

$$
W = \sum_{n=1}^{M} |\vec{F_{D_{n-1}}}| \cdot |\vec{r_n} - \vec{r_{n-1}}|,
$$
where $\vec{F_{D_{n-1}}}$ denotes the drag force at position $\vec{r_{n-1}}$. However, this can again be written to the form

$$
W = h \sum_{n=1}^{M} |\vec{F_{D_{n-1}}}| \cdot |\vec{v_{n-1}}|,
$$
which is the most convenient for our purposes. It must be noted that the choice of  index that is used for the force and speed respectively is non-trivial, as they both will vary during the interval $dr$. Here both the speed of the projectile and the force at point $\vec{r_{n-1}}$ is used. This is mainly a pragmatic choice, as this combination yields the lowest, and thus the most realistic energy loss. Thus we can finally express the energy of the projectile (which of course is equal to it's kinetic and potential energy) and the work done on the air at time $t_M$

$$
E(t_M) = h \sum_{n=1}^{M} |\vec{F_{D_{n-1}}}| \cdot |\vec{v_{n-1}}| + \frac{1}{2} m v_m^2 + mgy_m,
$$

where $v_m$ is the absolute velcoity of the projectile, and $y_m$ is the height of the projectile over the ground (zero potential reference point).

Now all that remains is to check if this expression for the energy indeed remains constant while it is in the air. 


In this notebook, the Runge-Kutta method is implemented without much difficulties. However, it is interesting to see whether or not it was necessary to use a 4th order differential equation solver, instead of using the most primitive ODE solver, Euler's method. The following code will use both Runge-Kutta and Euler's method, and compare the energy loss during the motion for both of them.


```python
def euler_step(t,w ,h, *P):
    """Simple implementation of Euler's method on vector form
    """
    return w + h * f(t,w)

def energy_conservation(h, theta, func):
    """ This function performs the same steps as "projectile_motion()" and "f()"
    but extracts the value for the drag, F_D, at every step. This enables the 
    calculation of work and storage of energy at each time step.
    Parameters:
        h           integration time step
        theta       initial shooting angle
        func        a function defining the integration step (euler or RK4)
    Returns:
        x_list      array of the trajectory x-coordinates
        y_list      array of the trajectory y-coordinates
        time_steps  array of the elapsed time at each step
        energy      array of the energy at each step
        work        the total work by the drag force
    """
    w = shoot(theta, v_0)
    
    x_list = np.zeros(1)
    y_list = np.zeros(1)
    work = 0
    
    # Initial total energy
    v = np.sqrt(w[2] ** 2 + w[3] ** 2)
    energy = np.array([0.5 * m * v ** 2])
    
    
    while w[1] >= 0:
        w = func(0, w, h,f,0)
        
        # Updating lists
        x_list = np.append(x_list, w[0])
        y_list = np.append(y_list, w[1])
        
        # Temporary parameters for checking if high altitude
        k = (1 - a * y_list[-2] / T_0)
        s = np.sqrt((w[2] - V[0]) ** 2 + (w[3] - V[1]) ** 2)
        if k > 0:
            F_D = B * k ** (alpha) * s ** 2
        else:
            F_D = 0
        
        # Add the amount of work done by drag this time interval
        work +=  v * F_D * h
        # Updating the speed in order to calculate the total energy
        v = np.sqrt(w[2] ** 2 + w[3] ** 2)
        # Use np.sum on the work_list to get the cumulative work
        energy = np.append(energy, 0.5 * m * v ** 2 + m * g * w[1] + work)
    
    n = len(x_list)
    time_steps = np.linspace(0, n * h, n)
    
    return x_list, y_list, time_steps, energy, work
 

# Setting initial values
theta = 45
h = 0.1
# Getting the results
x_list_euler, y_list_euler, time_steps_euler, energy_euler, work = energy_conservation(h, theta, euler_step)
x_list_RK4, y_list_RK4, time_steps_RK4, energy_RK4, work = energy_conservation(h, theta, RK4_step)

print("At firing angle %.1f degrees and with time step h = %.2f seconds:" %(theta, h))
print("Energy lost with RK4: %.3E = %.2f percent of total energy" %(energy_RK4[0] - energy_RK4[-1],
                                          (energy_RK4[0] - energy_RK4[-1]) / energy_RK4[0] * 100))
print("Energy lost with Euler's method: %.3E = %.2f percent of total energy" %(energy_euler[0] - energy_euler[-1],
                                        (energy_euler[0] - energy_euler[-1]) / energy_euler[0] * 100))

# Plotting the separate trajectories 
plt.plot(x_list_euler/1000, y_list_euler/1000, label ="Euler")
plt.plot(x_list_RK4/1000, y_list_RK4/1000, label ="RK4")
plt.title("Firing angle = " + str(theta) + "$\degree$")
plt.xlabel(r"$x$ [km]")
plt.ylabel(r"$y$ [km]")
plt.gca().set_ylim(bottom=0)
plt.legend()
plt.show()
```

From the figure above, it is clear that Euler's method work surprisingly well, even though it us just ODE-solver of first order. It produces nearly the same trajectory as Runge Kutta 4th method. The most likely explanation of the small difference between Runge Kutta 4th order and Euler's method is that the coupled set of ODEs we are solving are quite linear, in other words, the derivatives vary quite slowly. Thus Euler's method produce satisfactory results, even though it is only of first order.

It might seem like loosing around $10^3$ joules is a lot, but one needs to remember the enormous dimensions of the Paris gun, and that this is a small fraction of the original energy of the projectile.


```python

```

60615d454c2923296220c5f9baa336ef950a4e60
ipynb
Jupyter Notebook
Lectures/Lecture 13/Lecture13_ODE_part4.ipynb
astroarshn2000/PHYS305S20
18f4ebf0a51ba62fba34672cf76bd119d1db6f1e
MIT
2020-09-10T06:45:46.000Z
2020-10-20T13:50:11.000Z
Lectures/Lecture 13/Lecture13_ODE_part4.ipynb
astroarshn2000/PHYS305S20
18f4ebf0a51ba62fba34672cf76bd119d1db6f1e
MIT
Lectures/Lecture 13/Lecture13_ODE_part4.ipynb
astroarshn2000/PHYS305S20
18f4ebf0a51ba62fba34672cf76bd119d1db6f1e
MIT
Qwen/Qwen-72B
1. YES
2. YES