Instructions: Create a new notebook called ODEsExercises in your ODEs directory and solve the following problems inside it. Be sure to include the problem statements in a markdown cell above your solution. You don't need to put the "helper" code in the markdown cell, just implement the helper code in your code cell with your solution.
Preliminaries: At the top of your notebook, include a "Heading 1" cell with the title Ordinary Differential Equations Exercises. Then include the inline functions libraries by adding a code cell that invokes the %pylab inline magic and imports the needed packages.
For this question refer to the solutions for the nonlinear pendulum we got using scipy's odeint from the beginning of the ODEs Tour.
(a) Use the numerical solution to see what happens to the period of the pendulum as you change the initial angle of release. Plot $\theta$ vs. $t$ for 5 different initial angles, $\theta_0$ , spanning 0 to $\pi$. (i.e. 0, $\pi/4$ , $\pi/2$ , $3\pi/4$ , and $\pi$ ). Put them all on the same plot. What happens to the period of the pendulum as the initial angle increases from $\pi/4$ to $3\pi/4$? What happens to the trajectory when the initial angle is 0 or $\pi$ ?
In [1]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(b) Make plots of $\theta$ vs. $t$, $\omega$ vs. $t$, and $\omega$ vs. $\theta$ (the so-called "phase plot") for $\theta_0$ = 0.95$\pi$. Look carefully at their behavior. Explain what is going on in these figures.
In [2]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(c) For the $\theta_0$ = 0.95$\pi$ case, check that energy is conserved by plotting $E/m$ vs. $t$ (where $E/m =g\ell(1-\cos(\theta)) + \frac{1}{2}\ell^2\omega^2$ ).
Is the energy a constant? If not, then the computational errors are producing a non-physical loss (or gain) of energy. Review the documentation for odeint and change the error tolerance value to a tighter number.
In [3]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(d) Now look at what happens if we start the pendulum straight up, $\theta = \pi$. Go for at least 50 s. Be sure to describe what happens numerically according to the $\theta$ vs. $t$, $\omega$ vs. $t$ and $E/m$ vs. $t$ plots. What do you think should happen? What happens if you change the error tolerance on odeint? How does this compare to a nearly vertical pendulum:
$\theta_0 = \pi −0.000001$ ?
In [4]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(a) Use Euler's method (not odeint) to solve the following first-order differential equation with given initial condition $y(0)$ = 0.
Find the value of $y$ when $x$ = 1, $y(1)$, using two different step sizes: 0.25 then 0.02.
In [5]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(b) Plot the results from part (a) for the two different step sizes with the exact solution
$$ y = 0.25 e^{2x} - 0.5 x - 0.25 $$on one plot. Include a legend. Which step size gives a more accurate result?
In [6]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(c) Now solve the same problem using Euler-Cromer and the Midpoint method for the same pair of step sizes and plot the results along with the exact solution. Then also plot the difference between the approximate numerical solutions and the exact solution. Which combination of method and step-size gives the best result?
In [ ]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
We can use numerical solvers and widgets to explore the dynamics of a projectile that is subject to a drag force proportional to the square of its velocity - "quadratic drag".
The force on the projectile is given by
$$ \sum\vec{F} = m\vec{a}\\ = m\vec{g} - c v^2\hat{v}\\ = m\vec{g} - c v \vec{v}\\ = m\vec{g} - c \sqrt{v_x^2 + v_y^2}\vec{v} $$Where the second term is the drag force, in the direction opposite that of the velocity vector. Re-writing these in terms of the first-order differential equations for $v_x$ and $v_y$ gives
$$ m\frac{dv_x}{dt} = -c \sqrt{v_x^2 + v_y^2} v_x $$$$ m\frac{dv_y}{dt} = -mg -c \sqrt{v_x^2 + v_y^2} v_y $$taking the positive $y$ direction to be down. Notice that these are coupled equations. There are both $v_x$ and $v_y$ terms in each one. The drag coefficient $c$ depends on the shape and density of the projectile. We can re-write this system of equations in terms of the terminal velocity - the point at which the gravitational force balances the drag force. This happens when $mg = cv^2$. Let's denote the terminal velocity as $v_{ter} = \sqrt{mg/c}$ and re-write our equations in terms of $v_{ter}$.
$$ \frac{dv_x}{dt} = -g \frac{v_x}{v_{ter}^2} \sqrt{v_x^2 + v_y^2} $$$$ \frac{dv_y}{dt} = -g\left(1 + \frac{v_y}{v_{ter}^2} \sqrt{v_x^2 + v_y^2}\right) $$(a) Create a derivative function to solve this system of equations and use odeint to solve it for a projectile fired at an angle of 50 degrees from the horizontal at an initial speed of 30 m/s. Track the projectile's motion over the time interval of 0 to 200 seconds and let ($x_0$,$y_0$) = (0,0). You'll have to calculate ($v_{x0}$,$v_{y0}$) from the information provided. Let $v_{ter}$=44 m/s. Plot the trajectory of the projectile ($y$ vs. $x$) to see how far it can travel in the $x$ direction before landing.
In [8]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(b) Now create an interact object that uses sliders for the initial speed, terminal velocity, launch angle, and time. Have it plot two curves - one representing the solution with quadratic drag and the other in the absence of drag:
In [9]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(c) Use your interact object to find 3 pairs of initial angle and speed needed for each of the following projectiles to reach a target 80 m away. How do the speed values change if there is no air resistance? Put your results in a table for easier reading.
| Projectile | Terminal velocity (m/s) | Initial height (m) |
|---|---|---|
| Golfball | 44 | 0 |
| Baseball | 43 | 1.5 |
| Pumpkin | 107 | 1.5 |
In [11]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
The nonlinear pendulum is nice, but let’s make it more interesting by adding a damping and a driving force. This pendulum has much more interesting behaviors.
$$ \frac{d^2\theta}{dt^2} = \frac{-g}{\ell}\sin\theta - q\omega - f_d \sin(\Omega_d t) $$The damping force is governed by the factor $q$. The driving force is sinusoidal with an amplitude $f_d$ and angular frequency $\Omega_d$.
(a) Write a function and interact widget that solves the ODE for the nonlinear pendulum with damping and driving. The parameters, $q$, $f_d$, and $\Omega_d$ should be passed into the interact as sliders. Use initial values of $\theta_0$ =$\pi$/2, $\omega_0$=0 and a time span of 50 seconds.
Check that the solution looks as it did before for the case: $\ell$=0.5, $q$=0, $f_d$=0, $\Omega_d$=0. Play with the sliders to see how the solution depends on the parameters.
In [11]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell
(b) Make static plots of $\theta$ vs. $t$, $\omega$ vs. $t$, $E/m$ vs. $t$, and $\omega$ vs. $\theta$ for the ODE solution with
Include a legend with each plot. Do they look like you’d expect? Why or why not? Please explain what the results in each plot mean.
In [11]:
#Copy the exercise statement to a markdown cell in your notebook and then implement a solution in a code cell