**This is an old revision of the document!** ----
====== Math 445 Lab 11: Differential Equations ====== ===== Problem 1: dy/dt = sin(y+t) ===== Show a quiver plot in the //(t,y)// plane of the differential equation //dy/dt = sin(y+t)// over the region //0 ≤ t ≤ 2pi// and //-pi ≤ y ≤ pi//. Use Matlab's //ode45// function to compute solution curves //y(t)// for each of the following initial conditions, and superimposed the solution curves on the quiver plot in the indicated colors //y(0) = -1// in red, //y(0) = 0// in blue , //y(0) = 1// in black. Be sure to set the axes tight around the quiver plot, set the unit length of the axes to be equal, label the axes, and title the graph. ===== Problem 2: van der Pol oscillator ===== (a) The van der Pol oscillator is defined by the second-order differential equation <latex> d^2x/dt^2 - \mu (1-x^2) dx/dt + x = 0 </latex> This second-order equation can be converted to a system of two first-order equations by letting //y// be a 2-d vector $y = [y_1, y_2] = [x, ~dx/dt]$. Then $dy/dt = [dy_1/dt, ~dy_2/dt ] = [y_2, ~\mu (1-y_1^2) y_2 - y_1]$. Use this definition of the vector $dy/dt$ to produce a quiver plot of the van der Pol oscillator in the $y_1, y_2$ plane for $\mu = 1$, over the range $-3 \leq y_1 \leq 3$ and $-3 \leq y_2 \leq 3$. Using Matlab's //ode45//, compute solution curves $y(t)$ for the five different initial values $y(0) = [0.1~ 0], [0.2~0], [1~0], [2~0],$ and $[3~0]$, and superimpose them on the quiver plot in red. Hint: put the call to //ode45// inside a //for// loop that sets the value of $y_1$ to the five different initial values. (b) Put your work for part (a) in a script file, and run the script for a range of values of $\mu$ between -1 and 1. How does the character of the solutions of the equation change as $\mu$ passes through 0? Turn in plots for a negative, zero, and a positive value of $\mu$, whichever values you think best illustrate the change as $\mu$ passes through zero. Hint: To avoid the need to change the value of $\mu$ in both the script and the function file that defines $dy/dt = f(t,y)$, define $dy/dt = f(t,y)$ as an anonymous function within your script file, so that it inherits the $\mu$ of the script file. ===== Problem 3: projectile motion ===== It is easy to determine the trajectory of a projectile subject to gravity if you neglect air resistance. using elementary physics and calculus, you can show that a projectile will travel farthest across a flat surface if you shoot it at an angle $\theta = \pi/4$ above the horizontal. In this problem you will determine the optimal angle and the maximum range of a projectile subject to air resistance. Use the following equations of motion, developed during class $ dx/dt = v_x $ $ dy/dt = v_y $ $ dv_x/dt = -\mu v_x \sqrt{v_x^2 + v_y^2}$ $ dv_y/dt = -g - \mu v_x \sqrt{v_x^2 + v_y^2}$ The constant $g = 9.81 m/s^2$ is the acceleration due to gravity. The constant $\mu = 1/2 \rho_{air} C_D A/m$ in the air resistance term depends on physical characteristics of the projectile and the air. We will solve the problem for an iron cannonball 20cm in diameter. The following code will calculate $\mu$ correctly given appropriate physical constants <code> rho_air = 1.28; % kg/m^3, density of air rho_iron = 7870; % kg/m^3, density of iron 7.87 gm/cm^3 == 0.00787 kg/(0.01m)^3 = 7870 kg/m^3 C_D = 1; % drag coefficient for sphere g = 9.81; % m/s^2, acceleration due to gravity r = 0.10; % m, radius of cannonball A = pi*r^2; % cross-sectional area of cannonball m = rho_iron*4/3*pi*r^3; % mass of cannonball mu = rho_air*C_D*A/(2*m); % coefficient of nonlinear v|v| term, in mks units </code> Write a Matlab script to plot the trajectory of a cannonball subject to the above equations of motion and shot with an initial speed of 100 m/s at an angle $\theta$ above the horizontal. In the same plot, show the trajectory in the (x,y) plane of the cannonball with and without air resistance. Use your script to answer the following questions: (a) How far does the cannonball go in x with and without air resistance, if you shoot it at $\theta = \pi/4$? (b) What initial angle gives the farthest range with air resistance? (c) If $\theta = pi/4$, how much do you have to increase the initial velocity so the x distance is as large as it is without air resistance? Answering these questions with any precision will require some cleverness and a number of Matlab programming concepts you've learned through the semester.
