MATLAB LAB

matlab_lab2.pdf

Home >Mathematics homework help >MATLAB LAB

LAB #2

Escape Velocity

Goal: Determine the initial velocity an object is shot upward from the surface of the earth so as to never return; illustrate scaling variables to simplify di↵erential equations.

Required tools: dfield and pplane ; second order di↵erential equations; convert- ing a second order di↵erential equation with missing independent variable to a first order di↵erential equation; converting a second order di↵erential equation to a system of di↵erential equations.

Discussion

Newton’s Law of Gravitational Attraction tells us that in general an object of mass m which is x miles above the surface of the earth will feel a force of approximately

F = � 0.0061m

⇣ 1 +

4000

⌘2

pounds with the number “4000” an approximation to the radius of the earth in miles and “0.0061” the acceleration due to gravity measured in miles/sec2.

On the other hand, according to Newton’s 2nd Law

F = m d

2 x

2 .

Thus equating these expressions for F leads to the di↵erential equation

2 x

2 = �

0.0061 ⇣ 1 +

4000

⌘2 .

To study this equation using the Student Edition of Matlab, we need to change the units to avoid overly large intervals. Rather than measuring distance in miles, we will use multiples of the radius of the earth. This amounts to making the substitution

y = x

4000

(scaling the dependent variable) resulting in the equation

2 y

2 = �

0.0061

4000

(1 + y)2 . (A)

We can further simplify this equation by changing the units of time (scaling the independent variable). If we let

s =

r 0.0061

4000 t

then by the Chain Rule dy

= dy

r 0.0061

4000

and hence d

2 y

2 =

✓ dy

◆ =

r 0.0061

4000

2 y

2 =

r 0.0061

4000

✓ dy

◆

2 y

2 =

r 0.0061

4000

 d

✓ dy

◆ ds

�

2 y

2 =

0.0061

4000

2 y

Using the last expression, the di↵erential equation (A) becomes

2 y

2 = �

(1 + y)2 . (B)

Hence y is now viewed as a function of s, instead of t. Equation (B) is a second order di↵erential equation whose independent variable s

is missing. We can convert this special type of 2nd order di↵erential equation to a first order equation using a simple technique. Let

v = dy

This is simply the velocity of the projectile using s as the unit of time instead of t and y as the distance instead of x.

Since dy

= v, from the Chain Rule we obtain

2 y

2 =

= dv

= v dv

Thus, equation (B) can be written as

= � 1

(1 + y)2 . (C)

Equation (C) is a 1st order di↵erential equation which may be solved to express v as a function of y.

Assignment

(1) Enter equation (C) into dfield and plot a few solutions starting at y = 0 and various initial values for the velocity v using the dfield default ranges for v and y. You will need to enter it in the form

= � 1

v (1 + y)2 . (D)

(If you have trouble with dfield crashing, see the next exercise.)

(2) When the author did (D) using dfield7 with the Student Edition of Matlab, dfield crashed almost every time. The reason is that for this particular equation, the solution algorithm on which dfield is based requires many data points to produce an accurate picture. In fact, it can require so many as to exceed the maximum capacity of the Student Edition of Matlab. There is another program, pplane , which works much better for this equation. To use pplane , note that equation (B) is equivalent to the system

8 <

0 = v

0 = � 1

(1 + y)2

This is a system of di↵erential equations in 2 unknowns, y and v (s is the inde- pendent variable). You will study systems in detail later. For the moment, it is su�cient to know that pplane will plot one of the unknowns in the system as a function of the other. Specifically, start pplane (the command is “pplane ” or possibly “pplaneN” where N is the version of Matlab you are using) and then enter the above equations into the equation boxes of pplane . Use the ranges �1  y  15 and �2  v  3.

(3) Plot the solutions with initial values y = 0 and v = b where b is 0.5, 1.0, 1.5, and 2. You can either left click onto the appropriate points in the pplane “Display” window or use the “keyboard input” option in the “Solutions” pull down menu of the “Display” window. What is the maximum height of the object (in miles) in each case where there is a finite maximum ? (Recall that 4000y = x.) In the cases where there is no maximum, approximate the limiting value of v as s �! 1. (This is the terminal velocity). What is the terminal velocity in miles/sec ?

Hint : Since 4000y = x, we get dx

= 4000 dy

= 4000 v ds

(4) In Part (3), you should discover that the escape velocity lies between v = 1 and v = 1.5 (since object never reaches a finite maximum height, it never returns to earth). To estimate the escape velocity, plot the orbit that goes through the point (y, v) = (10, 0) and note the velocity v0 of this orbit corresponding to x = 0. This velocity will be somewhat below the escape velocity. Use the zoom feature to estimate v0 to within 2 places after the decimal. What is the value of v0 in miles/sec ?

(5) Solve the separable di↵erential equation (D) with the initial condition v(0) = v0 to obtain the following expression

y = 2

2 + (2 � v20) � 1 .

(a) According to this formula, for v0 = 1, what is the maximum value of y ? Explain your answer mathematically and demonstrate using a pplane plot.

(b) Use the above formula to prove that if v0 = p

2, then y tends to 1 as v approaches 0. The value v = 0 is then the terminal velocity. Check this using an appropriate pplane plot.

Hint : The escape velocity will be the smallest value of v0 for which y �! 1 as v tends to some number.

(6) The above calculations ignore the e↵ect of air resistance on the object. We assume that resistance is proportional to velocity and decreases with increasing height. Explicitly, we assume the frictional force to be

Ff = � kv

(1 + y)2

where k is a physical constant (the coe�cient of friction) and we will initially take k = 0.01. We will also assume that our object has mass m = 1 in which case our di↵erential equation becomes

2 y

2 = �

(1 + y)2 �

0.01v

(1 + y)2 = �

(1 + 0.01v)

(1 + y)2 . (E)

Rewrite (E) as a system of equations as in Part (2). Use pplane to estimate the escape velocity as in Part (2) to find how much air resistance changes the calcu- lations. Next, solve the equation analytically. You will not be able to explicitly solve for v. Instead, you should get

100v � 104 ln(1 + 0.01v) � 1

(1 + y) = C

Show that C = 100v0 � 104 ln(1 + 0.01v0) � 1 and then explain why the object escapes if and only if 100v0 � 104 ln(1 + 0.01v0) > 1. Exit pplane and then use Matlab to plot the function 100v�104 ln(1+0.01v)�1 to find the escape velocity. Make sure that your answer is accurate to within two digits after the decimal. This may require doing several plots over increasingly smaller intervals. The appropriate sequence of commands is:

>> grid on

>> hold on

>> fplot(’100⇤x�10^4⇤log(1 + 0.01⇤x)�1’,[xmin,xmax,ymin,ymax]) where xmin, xmax and ymin, ymax are the ranges on the horizontal and vertical axes, respectively. Note that the function to be plotted must be enclosed in single quotes. The command “grid on” causes the axes to be visible and “hold on” prevents them from being erased as subsequent graphs are drawn.

(7) Estimate the escape velocity as in Part (6) if k = 0.1

seed