Math 449 - Homework 11 - Solutions
                   
                   
                   

1) Exercise 4, section 9.1, page 463.

(a) Equation: y' = exp(-2*t) - 2*y
    
    Solution: y(t) = (C + t)*exp(-2*t) 
    
Direct verification: y'(t) = exp(-2*t) - 2*(C + t)*exp(-2t)
 
                           = exp(-2*t) - 2*y(t)
                           
(b) Find a Lipschitz constant L for the rectangle R = [0, 3]x[0, 5]

By theorem 9.1, it is sufficient to find a constant L that bounds

from above the absolute value of the partial derivative of f(t, y) in y. 

This partial derivative is f_y(t,y) = -2. Therefore, we can take L=2.                              
    
    
************************************************************************** 

2) Exercise 6, section 9.1, page 463.

Construct a graph of the slope field m(i,j) = f(t_i, y_j) over the rectangle

[0, 4]x[0, 4] for the differential equation y' = -t/y and draw the solution

curves y(t) = (C - t^2)^(1/2) for C= 1, 2, 4, 9. Note that it will

be necessary to restrict the range of values of t in each case so that

the solution exists over that range. The interval of t is [0, sqrt(C)] in 

each case. 


I will use the script shown on top of page 463 (adapted for this equation):

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[t, y] = meshgrid(.5:.5:3.5, .5:.5:3.5);
dt = ones(7,7);
dy = -t./y;
quiver(t, y, dt, dy);
hold on
x1 = 0:.01:1;
x2 = 0:.01:sqrt(2);
x3 = 0:.01:2;
x4 = 0:.01:3;
z1 = (1 - x1.^2).^(1/2);
z2 = (2 - x2.^2).^(1/2);
z3 = (4 - x3.^2).^(1/2);
z4 = (9 - x4.^2).^(1/2);
plot(x1, z1)
plot(x2, z2)
plot(x3, z3)
plot(x4, z4)
grid
axis equal
hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

See the graphs on the separate page.

************************************************************************** 

3) Exercise 7, section 9.2, page 473.

Show that when Euler's method is used to solve the I.V.P.

           y' = f(t) over [a, b] with y(a) = y_0 = 0
           
the result is the Riemann sum that approximates the definitive integral

of f(t) taken over the interval [a, b].


The recursive formula for the Euler iteration is 

y_(k+1) = y_k + f(t_k, y_k)*h, which in this case reduces to

y_(k+1) = y_k + f(t_k)*h. The total sum is

y_M = y_0 + f(t_0)*h + f(t_1)*h + ... + f(t_(M-1))*h,

which is the Riemann sum that approximates y(b) = Int_[a,b] f(t) dt.
           


3') This is problem 7, of Algorithm and Programs, 9.2, page 474, which

I did accidentally, instead of the previous one. 

Exponential polulation growth. The population of a certain species grows

at a rate that is proportional to the current population and obeys the

I.V.P.

              y' = 0.02*y  over  [0, 5]  with  y(0) = 5000.
              
(a) Apply the formula y_k = y_0*(1 + h*R)^k, t_k = kh, k = 0, 1, ..., M

to find Euler's approximation to y(5) using the step sizes h= 1, 1/12, 1/360.

M=5, 5*12, 5*360.
  
Case 1: M = 5,    y(5) approx. 500*(1 + 0.02)^5            = 552.0404

Case 2: M = 60,   y(5) approx. 500*(1 + (1/12)*0.02)^60    = 552.5395

Case 3: M = 1800, y(5) approx. 500*(1 + (1/360)*0.02)^1800 = 552.5839


(b) The limit is y(5) = 500*exp(0.02*5) = 552.5855
************************************************************************** 

4) Algorithms and Programs 8, section 9.2, page 474.

A skydiver jumps from a plane, and up to the moment he opens the parachute

the air resistance is proportional to v^(3/2). (v represents velocity).

Assume that the time interval is [0, 6] and that the differential equation

for the downward direction is 

              v' = 32 - 0.032*v^(3/2) over [0, 6] with v(0) = 0.
              
Use Euler's method with h = 0.05 and estimate v(6). (Note: h = 0.05 implies

M = 6/0.05 = 120.)


The following program implements Euler's method:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function E=euler(f,a,b,ya,M)
%Input  - f is the function entered as a string 'f'
%       - a and b are the left and right endpoints
%       - ya is the initial condition y(a)
%       - M is the number of steps
%Output - E=[T' Y'], where T is the vector of abscissas and
%         Y is the vector of ordinates.
h=(b-a)/M;
T=a:h:b;
Y=zeros(1,M+1);
Y(1)=ya;
for j=1:M
    Y(j+1)=Y(j)+h*feval(f,T(j),Y(j));
end
E=[T' Y'];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I denoted the slope field function g.m:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=g(t,x)
y=32-0.032*x.^(3/2);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

With the above saved as .m files, we now call the command

E=euler('g',0,6,0,120)

This gives the matrix whose first column is the time values and the second

column the corresponding distance moved downward from the place of jump.

The velocity at time 6 is 92.4979 (units were not specified in the problem,

but it should be seconds and feet since the acceleration of gravity in 

imperial units is 32 feet per square second.)
************************************************************************** 

5) Exercise 6, section 9.3, page 480.

Show that when Heun's method is used to solve the I.V.P. y' = f(t) over

the interval [a, b] with y(a) = y_0 = 0 the result is the trapezoidal

rule approximation for the definite integral of f(t) taken over [a, b].


The Heun's iteration has the form

p(k+1) = y(k) + h*f(t(k), y(k))

y(k+1) = y(k) + (h/2)*[f(t(k),y(k)) + f(t(k+1), p(k+1))]

If f does not depend on y, this system reduces to

p(k+1) = y(k) + h*f(t(k))

y(k+1) = y(k) + (h/2)*[f(t(k)) + f(t(k+1))]

The first equation becomes irrelevant and the second, summed over k, gives

y(M) = y(0) + (h/2)*[f(t(0)) + f(t(1))] + ... + (h/2)*[f(t(M-1)) + f(t(M))]

Taking y(0) = 0, the result is the trapezoidal rule approximation for 

y(M) = Int_[a, b] f(t) dt.
************************************************************************** 

6) Algorithms and Programs 6, section 9.3, page 482.

Consider a projectile that is fired straight up and falls straight down. If

air resistance is proportional to the velocity, the I.V.P. for the velocity

v(t) is  v' = -32 - (K/M)*v  with  v(0) = v_0 where v_0 is the initial 

velocity, M is the mass, and K the coefficient of air resistance. Suppose

that v_0 = 160 ft/sec and K/M = 0.1. Use Heun's method with h = 0.5 to solve

             v' = -32 - 0.1*v  over  [0, 30]  with v(0) = 160.
             
Graph the solution and the exact solution v(t) = 480*exp(-t/10) - 320 on the

same coordinate system.

The following program implements Heun's method:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function H=heun(f,a,b,ya,M)
%Input  - f is the function entered as a string 'f'
%       - a and b are the left and right endpoints
%       - ya is thte initial condition y(a)
%Output - H=[T' Y'], where T is the vector of abscissas and
%         Y is the vector of ordinates
h=(b-a)/M;
T=a:h:b;
Y=zeros(1,M+1);
Y(1)=ya;
for j=1:M
    k1=feval(f,T(j),Y(j));
    k2=feval(f,T(j+1),Y(j)+h*k1);
end
H=[T' Y'];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The slope field is written on the function file g.m:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=g(t,x)
y=-32-0.1*x;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

I will use h = 0.05, so M=30/0.05=600. The command to run the function 

heun is: H=heun('g', 0,30,160, 600);

To plot, write T=H(:,1) and V=H(:,2) then plot(T,V)

The graph is shown on a separate page.
**************************************************************************