=================================================== Math 451 Laboratory 4: Differential Equations Date: 2013-04-05 Due: 2013-04-12 =================================================== Differential equations are one of the most useful conceptual tools of applied mathematics. But solving them requires some practice. Goals ---------------------------- - Learn how to solve a differential equation with Matlab - Learn to nth order equations as systems - Learn to plot time series and phase portraits - Solve Newton's equations for a modified 3-body problem - Make javascript animations of differential equation - Solve Newton's planet equations with python Note: The differential equation solver in Octave (lsode) is not compatible with Matlab, so you will have to do some creative re-writing as you go if you would like to complete this lab using Octave. =================================================== A simple example of an initial value problem ---------------------------------------------- The first-order equation dn/dt = r n (1 -n ) - k (1 - cos(2*pi*t) ) is what's called a Ricatti equation. The variable n(t) represents the density of oysters relative to their natural carrying capacity and the parameter r represents the maximum per-capital growth rate of the population. The oysters compete for space on the ocean bottom, and they are harvested at different rates over the course of the year. We would like to know if harvesting rates corresponding to k = 0.2 are too intense and will collapse the stock. To do this, we'd like to explore the solutions of this equation, but unfortunately, this equation does not have closed-form solutions. So instead, we will solve it numerically. To solve this equation numerically, we need to specify the equation's right-hand side, an initial condition, and a sequence of time-points where we would like to calculation the value of the solution curve. (1-a) Create a Matlab function file called "oyster_ode.m" This function takes the current time and state as arguments and returns a vector representing the direction and rate of change in the next time-step. function dn = oyster_ode(t,n) k = 0.2; r = 1.0; dn = [r*n(1)*(1 - n(1))-k*(1-cos(2*pi*t))]; return (1-b) Create a Matlab script called 'oyster_simulator.m' using the code below. ic = [ 1. ]; tspan = linspace(0,8,1000); [t,x] = ode45( 'oyster_ode', tspan, ic); figure(1); plot(t,x,'o-'); xlabel('Time'); ylabel('Oyster density'); ylim([-1,3]); This file specifies and solves the differential equation. Here, ic sets the initial oyster population size and tspan is the initial and final times at which we would like to know solution. ode45 is a numerical method for solving differential equations. Specifically, it is a 5th order Runge-Kutta method with adaptive time stepping. It returns a vector t of times, and a matrix x of the solution at those times. (1-c) Run your script 'oyster_simulator.m'. (1-d) Does the stock crash and go extinct? (1-e) To show how what happens depends on the initial population size, plot solutions for 20 initial conditions ranging from 0 to 3 on the same plot. Save your figure to be handed in. You may find the "hold on" command useful. Solutions of systems -------------------------------------------------------- Numerical solvers like ode45 can only solve 1st-order equations. So whenever we need to solve a 2nd order equation, we transform it into a system of first order equations and solve those. Consider the damped pendulum equation u" = - c u' - sin( u ) with c = 0.1. This is a second order nonlinear equation. (2-a) Use the substitution y(1) = u, y(2) = u' to rewrite this equation as a system of 2 first-order equations. (2-b) Write a new function file called "pendulum_vfield.m" in the same style as "oyster_ode.m" that takes a time and a position y as input arguments and returns a vector representing the current state velocity AS A COLUMN VECTOR. (2-c) Write a new script named "pendulum_simulator.m" that plots time-series for u(t) and u'(t) starting from an initial position u(0) = pi/2 and initial velocity u'(0) = 0 over 25 seconds. Include labels and a legend. Save this figure to hand in. Phase portraits are plots of the states of 2-dimensional systems relative to each other. For some simple examples, see http://tutorial.math.lamar.edu/Classes/DE/PhasePlane.aspx and http://en.wikipedia.org/wiki/Phase_portrait (2-d) Construct solutions to the pendulum equation with u'(0) is in linspace(-5,5,15) and u(0) = 0, and plot these solutions to make a phase portrait. Label your axes and save this figure to hand in. (2-e) Why are the spirals repeating? (optional) What are the equilibrium points and where are they in the phase portrait? What are the basins of attraction of attracting equilibria? Javascript and Differential Equations in Browsers ---------------------------------------------------- One downside of the images we've produced so far is that they are static, and don't communicate any of the dynamic nature of differential equations. But our web browsers with javascript and HTML5 are now very powerful tools for doing this. For example, see http://hint.fm/wind/ Let's do a simpler version of this (the original version of this is at https://gist.github.com/1320178). (3-a) Download the "lorenz_equation.html" to your Desktop, open it, and read it. The file contains a combination of html markup, style-sheet specifications, and javascript code. * The first part in the "style" tags is the style sheet stuff * The next part is the script to solve our differential equation. Even though we haven't used javascript in this class before, you should be able to recognize variable names and equations to get an idea for how the code works. (3-b) Reconstruct the specification of the differential equations this code is solving using Euler's method. (3-c) What is the step size being used in Euler's method? Now, test the code in a web browser to see what it does. (3-d) You may notice that the dots change size. What does this mean? (3-e) How many particles are being used in the paths? Now, we'll make some simple changes to the code to change the animations we get out. (3-f) How does the behavior of the animation change as you use large step sizes? Be specific. (Don't increase the step size too fast, or it will be hard to tell.) (3-g) Modify the code so that we can see the cumulative shape of the attractor. Describe your changes, and include a screen-shot of the result. Newton's equations for planetary motion in python ---------------------------------------------------- The original killer application of differential equations was explaining planetary motion. Unfortunately, these equations are non-linear and at least 4-dimensional. This is too hard for us to solve in our first differential equations class, but they are not too hard to solve numerically! The main concept is that forces from gravity on a body are proportional to the masses of both bodies and inversely proportional to the distances between the bodies. If r is the radial position vector of a planet from it's sun, then Newton's equation is 2 d r G m r ---- = ------- 2 3 dt |r| (4-a) Rewrite this equation as a system of 4 1st-order differential equations where r = (x,y). We will solve these equations using a python script. (4-b) Download "planetary_motion.py". This script predicts the motional of a small satellite in a system with a star (or starts) that have fixed positions. (4-c) Run the script "planetary_motion.py" in the terminal using the command python planetary_motion.py One figure you get should be the path of a planet around 1 sun. The other is the time-series for the x and y positions and the x and y velocities. Does the picture agree with your and Johannes Kepler's expectations? (4-d) Add a second sun to the system, with the same mass as the first sun, but at position [-0.5,0.]. Keep everything else the same. (4-e) Run your new script and save the orbital path figure to hand in. (Extra credit) Port this script to Matlab.