Math 449 - Homework 12 - Solutions
                   
                   
                   

1) Algorithms and Programs 9, section 9.5, page 504.

By solving the appropriate initial value problem, make a table of values

of the function f(t) given by the following integral:

           f(t) = 1/2 + (1/sqrt(2*pi))Int_[0,x] exp(-t^2/2) dt
           
for x in [0, 3]. Use Runge-Kutta method of order N=4 with h = 0.1.


The initial value problem is y' = (1/sqrt(2*pi))*exp(-t^2/2), y(0) = 1/2.

The problem that implements Runge-Kutta method of order 4 is:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function R=rk4(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 - R=[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=h*feval(f,T(j),Y(j));
    k2=h*feval(f,T(j)+h/2,Y(j)+k1/2);
    k3=h*feval(f,T(j)+h/2,Y(j)+k2/2);
    k4=h*feval(f,T(j)+h,Y(j)+k3);
    Y(j+1)=Y(j)+(k1+2*k2+2*k3+k4)/6;
end
R=[T' Y'];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
           
We save the function to be integrated in the file named g.m:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function y=g(t,x)
y=(1/sqrt(2*pi))*exp(-t.^2/2); 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

We now choose a value for x, say x=3, and call the rk4 program as follows:

R=rk4('g',0,x,0.5,x/0.1);

This produces the table:


    x  |        f(x)
   ____|__________________        
   0.0 |  0.50000000000000
   0.1 |  0.53982784140429
   0.2 |  0.57925971748903
   0.3 |  0.61791143376648
   0.4 |  0.65542175615701
   0.5 |  0.69146247810558
   0.6 |  0.72574690060114
   0.7 |  0.75803636684998
   0.8 |  0.78814462042910
   0.9 |  0.81593989288217
   1.0 |  0.84134476288708
   1.1 |  0.86433395395856
   1.2 |  0.88493034240682
   1.3 |  0.90319952554994
   1.4 |  0.91924334833413
   1.5 |  0.93319280378551
   1.6 |  0.94520071100440
   1.7 |  0.95543453784294
   1.8 |  0.96406967969267
   1.9 |  0.97128343753301
   2.0 |  0.97724986429239
   2.1 |  0.98213557490590
   2.2 |  0.98609654749230
   2.3 |  0.98927588479114
   2.4 |  0.99180245891943
   2.5 |  0.99379032972536
   2.6 |  0.99533880736310
   2.7 |  0.99653302200441
   2.8 |  0.99744486594447
   2.9 |  0.99813418345723
   3.0 |  0.99865009919929

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

In the following problems we use program 9.9 which implements Runge-Kuta

method of order 4 for systems of differential equations of first order.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [T,Z]=rks4(F,a,b,Za,M)
%Input  - F is the system input as a string 'F'
%       - a and b are the endpoints of the interval
%       - Za=[x1(a) ... xn(a)] are the initial conditions
%       - M is the number of steps
%Output - T is vector of steps
%       - Z=[x1(t) ... xn(t)], where xk(t) is the approximation
%         to the kth dependent variable.
h=(b-a)/M;
T=a:h:b;
Z=zeros(M+1,length(Za));
Z(1,:)=Za;
for j=1:M
    k1=h*feval(F,T(j),Z(j,:));
    k2=h*feval(F,T(j)+h/2,Z(j,:)+k1/2);
    k3=h*feval(F,T(j)+h/2,Z(j,:)+k2/2);
    k4=h*feval(F,T(j)+h,Z(j,:)+k3);
    Z(j+1,:)=Z(j,:)+(k1+2*k2+2*k3+k4)/6;
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

2) Algorithms and Programs 12, section 9.7, page 526.

At time t, a pendulum makes an angle x(t) with the vertical axis. Assuming

that there is no friction, the equation of motion is

                     m*l*x''(t) = -m*g*sin(x(t))

where m is the mass and l is the length of the string. Use the Runge-Kutta 

method to solve the differential equation over the interval [0,2] using

M=40 steps and h=0.05 if g = 32 ft/sec^2 and

(b) l = 0.8 ft and x(0) = 0.3 and x'(0) = 0.


We first write this as a system of first order equations. Define:

x_1 = x  and X_2 = x' Then, using g/l = 10 sec^(-2),

                        x_1' = x_2
                        
                        x_2' = -10*sin(x_1)

with initial conditions: x_1(0) = 0.3, x_2(0) = 0.

Now apply rks4 (above) to the function:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U=g(t,Z)
U=[Z(2), -10*sin(Z(1))];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[T,Z]=rks4('g',0,2,[0.3 0],40);
plot(T,Z(:,1))

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

3) Algorithms and Programs 13, section 9.7, page 526.

We wish to solve the system

           x' = A*x - B*x*y
           
           y' = C*x*y - D*y

with parameter values A = 2, B = 0.02, C = 0.0002, D = 0.8.

The time interval is [0,5], M=50 steps and h=0.2. The initial conditions

are:

(a) x(0) = 3000 rabbits and y(0) = 120 foxes


The function file is g.m:

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U=g(t,Z)
A=2;
B=0.02;
C=0.0002;
D=0.8;
U=[A*Z(1) - B*Z(1)*Z(2), C*Z(1)*Z(2) - D*Z(2)];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

To obtain the graph, use

[T,Z]=rks4('g',0,5,[3000 120],50);
plot(Z(:,1),Z(:,2))

Notice how the numbers of both rabbits and foxes vary periodically.

When there are few foxes and rabbits, the number of rabbits increases

rapidly and the number of foxes slowly accelerates. As the fox population

increases beyond a certain point the number of rabbits begin to decrease 

until it reaches a point that cannot sustain the large fox population. 

The latter then decreases relatively fast until the two are back at the 

initial numbers.

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

4) Algorithms and Programs 17, section 9.7, page 528.

Solve x' = 1 - y,  y' = x^2 - y^2  over [0,5] using 0.1 (M = 50). 

Use the initial conditions (x(0), y(0)) = (2,0), (0,2), (0,0), (-1.2,0), 

and (-0.5,2).


Once again, we use the above rks4. The function is: 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function U=g(t,Z)
U=[1-Z(2), Z(1)^2-Z(2)^2];
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The graphs are shown on a separate page. 

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