The MATLAB Notebook v1.5



MATLAB 5 ODE Tutorial Using Baseball

C. Henry Edwards

The University of Georgia

hedwards@math.uga.edu

math.uga.edu/~hedwards

We illustrate some of the MATLAB 5 ODE suite's new facilities using the example of a batted baseball. With air resistance proportional to the square of its velocity v, the ball's horizontal and vertical coordinate functions [pic] and [pic] satisfy the acceleration equations

[pic] (1)

where c is a proportionality constant and g = 32 ft/sec2 is the downward acceleration of gravity. For computational purposes we convert (1) to the equivalent first-order system

[pic] (2)

where [pic], [pic] denote the ball's coordinate functions and [pic], [pic] its velocity components.

WARNING: The MATLAB functions that is typed below must be extracted and saved as separate m-files before the other commands in this notebook can be executed.

The following MATLAB function encodes the equations in (2), and is a simple example of an "ODE file".

type bbode0

function yp = bbode0(t,y)

% yp = bbode0(t,y)

% defines the first order baseball system

% with air resistance proportional to v-squared.

g = 32;

c = 0; % No air resistance yet

yp = y;

v = sqrt(y(3)*y(3) + y(4)*y(4));

yp(1) = y(3);

yp(2) = y(4);

yp(3) = -c*v*y(3);

yp(4) = -c*v*y(4) - g;

This function bbode0 takes as input the time t and the column vector

y = [y1 y2 y3 y4]T and outputs the column vector yp of derivatives.

The Old Fashioned Way

First let's review some of the standard textbook numerical ODE methods. The following function implements the Euler method with n equal steps of length [pic] from the initial time t with initial value y to time tfinal. It outputs a column vector T of successive times and a matrix Y whose rows are the corresponding approximate y-vectors.

type euler

function [T,Y] = euler(t,y,tfinal,n)

% The classical Euler method

h = (tfinal - t)/n; % step size

T = t; % initial t

Y = y'; % initial y-vector

for i = 1 : n % begin loop

k = f(t, y); % slope vector

t = t + h; % new t

y = y + h*k; % new y

T = [T;t]; % update t-column

Y = [Y;y']; % update Y-matrix

end % end loop

Here f denotes the ODE file to be used. The Euler method serves to illustrate the basic motif

[pic] (3)

of numerical ODE approximation, but is not sufficiently accurate for practical use.

The classical Runge-Kutta method uses an average of 4 carefully selected intermediate slopes for the "average slope" in (3). The following MATLAB Runge-Kutta function uses the same notation as in the Euler method above.

type rkn

function [T,Y] = rkn(t,y,tfinal,n)

% The classical Runge-Kutta method

h = (tfinal - t)/n; % step size

T = t; % initial t

Y = y'; % initial y-vector

for i = 1 : n % begin loop

k1 = f(t, y); % first slope

k2 = f(t+h/2, y+h*k1/2); % second slope

k3 = f(t+h/2, y+h*k2/2); % third slope

k4 = f(t+h , y+h*k3 ); % fourth slope

k = (k1+2*k2+2*k3+k4)/6; % average slope

t = t + h; % new t

y = y + h*k; % new y

T = [T;t]; % update t-column

Y = [Y;y']; % update Y-matrix

end % end loop

To illustrate the standard textbook approach, consider a batted baseball with initial position (0, 0), initial velocity v0 = 160 ft/sec, and initial inclination angle

( = 30(. According to the elementary exact solution it should remain in the air exactly [pic] seconds, and hit the ground at a distance of

[pic] feet away. The initial vector is then

[0; 0; 80*sqrt(3); 80] and if we make a copy of the file bbode0.m named f.m for rkn to use,

type f

function yp = f(t,y)

% yp = bbode0(t,y)

% defines the first order baseball system

% with air resistance proportional to v-squared.

g = 32;

c = 0; % No air resistance yet

yp = y;

v = sqrt(y(3)*y(3) + y(4)*y(4));

yp(1) = y(3);

yp(2) = y(4);

yp(3) = -c*v*y(3);

yp(4) = -c*v*y(4) - g;

then the command

[T,Y] = rkn(0,[0;0;80*sqrt(3);80],5,50);

yields a matrix Y of size

size(Y)

ans =

51 4

whose 51st row

Y(51,:)

ans =

692.8203 0.0000 138.5641 -80.0000

shows good accuracy in this simple problem. However, fixed stepsize methods like Runge-Kutta are not sufficiently accurate for many applications.

The MATLAB Way

Modern numerical ODE methods uses sophisticated variable-stepsize methods in which the stepsize is optimized as the solution progresses. The MATLAB ode45 solver implements a Dormand-Prince 4th-5th order method and has the basic syntax ode45('odefile', tspan, y0) where odefile (like bbode)defines the system being solved, tspan = [a b] is the desired time interval, and y(a) = y is the initial vector.

For our baseball problem with no air resistance we can enter

tspan = [0 5];

v0 = 160; alpha = pi/6;

y0 = [0; 0; v0*cos(alpha); v0*sin(alpha)]; % [x0; y0; x0'; y0']

[T, Y] = ode45('bbode0', tspan, y0);

How many steps were taken?

size(Y)

ans =

73 4

So the ball's position after 5 seconds is

Y(73,[1 2])

ans =

692.8203 0.0000

Thus, as expected, it hits the ground just under under 693 feet away. How accurate is the result?

ans - [400*sqrt(3), 0]

ans =

1.0e-012 *

-0.1137 0.0426

To about 10-13 foot! Let's plot this massive home run.

plot( Y(:,1), Y(:,2) ), axis([ 0 693 0 200

])

Coding Initial Conditions in the ODE file

It's a bit awkward to key in the initial velocity components, so let's encode the desired t-interval and the initial conditions.

type bbode1

function [out1,out2,out3] = bbode1(t,y,flag)

% [out1,out2,out3] = bbode1(t,y,flag)

% returns initial conditions

if nargin < 3 | isempty(flag)

% Return dy/dt = f(t,y)

g = 32;

c = 0;

yp = y;

v = sqrt( y(3)*y(3) + y(4)*y(4));

yp(1) = y(3);

yp(2) = y(4);

yp(3) = -c*v*y(3);

yp(4) = -c*v*y(4) - g;

out1 = yp;

elseif strcmp(flag,'init')

% Return [tspan, y0, options]

out1 = [0; 5]; % tspan

v0 = 160; alpha = pi/6; % initially

out2 = [0; 0; v0*cos(alpha); v0*sin(alpha)];

out3 = [];

end

Now we can simply say

[T,Y] = ode45('bbode1');

size(Y)

ans =

73 4

Y(73,[1 2])

ans =

692.8203 0.0000

Same results as before.

Including a Parameter in the ODE file

We can also include the air resistance proportionality constant c as a parameter in our ODE file that can be specified when we call ode45.

type bbode2

function [out1,out2,out3] = bbode2(t,y,flag,c)

% [out1,out2,out3] = bbode2(t,y,flag,c)

% returns initial conditions if flag = 'init'

if nargin < 4 | isempty(c)

c = 0; % default is no air resistance

end

if nargin < 3 | isempty(flag)

% Return dy/dt = f(t,y)

g = 32;

yp = y;

v = sqrt( y(3)*y(3) + y(4)*y(4));

yp(1) = y(3);

yp(2) = y(4);

yp(3) = -c*v*y(3);

yp(4) = -c*v*y(4) - g;

out1 = yp;

elseif strcmp(flag,'init')

% Return [tspan, y0, options]

v0 = 160; % default initial velocity

alpha = pi/6; % default init inclination

out1 = [0; 5]; % tspan

out2 = [0; 0; v0*cos(alpha); v0*sin(alpha)];

out3 = []; % options

end

[T, Y] = ode45('bbode2');

size(Y)

ans =

73 4

Y(73, [1 2])

ans =

692.8203 0.0000

The familiar result. But the elementary exact solution says we'd get the same range 1600 sin ( cos ( with ( = 60( as with ( = 30(, but with the ball hitting the ground after 10 sin ( = 5(3 (rather than 5) seconds. Let's check this out.

[T,Y] = ode45('bbode2',[0;5*sqrt(3)],[0;0;160*cos(pi/3);160*sin(pi/3)]);

size(Y)

ans =

77 4

Y(77, [1 2])

ans =

692.8203 0.0000

Yep!

Taking Air Resistance into Account.

Now let's suppose velocity-squared air resistance coefficient c = 0.0025. The ball shouldn't be in the air so long, so let's try 4 seconds.

[T, Y] = ode45('bbode2',[0;4],[],[],0.0025);

size(Y)

ans =

73 4

Y(73, [1 2])

ans =

339.6715 0.9087

Good guess! After 4 seconds the ball is just under a foot off the ground.

[T,Y] = ode45('bbode2',[0;4.1],[],[],0.0025);

Y(73, [1 2])

ans =

345.0965 -4.8392

So interpolation suggests the ball has traveled a bit more than 340 feet in just over 4 seconds. The air resistance has converted a massive home run into a routine fly ball (if hit to straightaway centerfield).

Event Location

Most real-world ODE problems call not merely for the approximation of a given ODE system with specified initial conditions, but for the location of some specific "event" associated with the solution. For instance, suppose that, instead of relying upon trial-and-error or interpolation, we want ode45 to automatically tell us when and where the ball hits the ground – and also, for good measure – when and where it reaches the apex of its trajectory. This is a matter of locating events – when y = 0 and when vy = 0. The following ODE file illustrates the definition of events to be located. In this version we also add v0 and alpha as parameters, as well as the exponent k in the air resistance model R = cvk.)

type bbode3

function [out1,out2,out3] = bbode3(t,y,flag,v0,alpha,c,k)

% [out1,out2,out3] = bbode3(t,y,flag,c,k)

% returns apex and impact events if flag = 'events'

if nargin < 7

k = 2; % default is v-squared

end

if nargin < 6 | isempty(c)

c = 0; % default is no air resistance

end

if nargin < 5 | (isempty(v0)&isempty(alpha))

v0 = 160;

alpha = pi/6;

end

if nargin < 3 | isempty(flag)

% Return dy/dt = f(t,y)

g = 32;

yp = y;

v = sqrt( y(3)*y(3) + y(4)*y(4));

r = c*v^(k-1); %resistance/v

yp(1) = y(3);

yp(2) = y(4);

yp(3) = -r*y(3);

yp(4) = -r*y(4) - g;

out1 = yp;

else

switch(flag)

case 'init' % Return [tspan, y0, options]

% v0 = 160; % default initial velocity

% alpha = pi/6; % default init inclination

out1 = [0; 10]; % max tspan

out2 = [0; 0; v0*cos(alpha); v0*sin(alpha)];

out3 = odeset('Events','on','RelTol',1e-4);

% options

case 'events'

out1 = [y(4); y(2)]; % 1st zero event is vy = 0

% 2nd zero event is y = 0

out2 = [0; 1]; % Don't stop when vy = 0

% DO stop when y = 0

out3 = [-1;-1]; % Actual event when either is

% decreasing through zero

otherwise

error(['Unknown flag ''' flag '''.']);

end

end

Let's check it out, first with the default c = 0 (no air resistance).

[T,Y,te,ye,ie] = ode45('bbode3');

In this extended output,

te

te =

2.5000

5.0000

specifies the times at which the apex and impact events occur, and

ye

ye =

346.4102 100.0000 138.5641 0.0000

692.8203 0.0000 138.5641 -80.0000

gives the y-vectors at these events.

Thus:

Thus the first event occurred at time t = 2.5000 when the ball was 100 ft high.

The second event occurred at time t = 5.000 when the ball had traveled 692.8203 feet.

Now let's see what happens with air resistance coefficient c = 0.0025 (but default values for the other parameters)

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],0.0025);

te

te =

1.8532

4.0161

ye

ye =

194.5579 65.6439 82.4561 0.0000

340.5509 0.0000 54.6148 -56.7772

Thus:

Thus the first event occurred at time t = 1.8532 when the ball was only 65.64 feet high.

The second event occurred at time t = 4.0161 when the ball had traveled 340.55 feet.

We might note that the ball hits the ground with velocity

sqrt( ye(2,3)^2 + ye(2,4)^2 )

ans =

78.7811

having lost to air resistance over half its original v0 = 160 ft/sec. It impacts at

angle

180*atan(-ye(2,4)/ye(2,3))/pi

ans =

46.1123

Thus at 46 degrees, steeper than its initial inclination of 30 degrees. Every baseball fan has observed empirically these aspects of the trajectory of a fly ball.

Shooting (i.e., batting) Problems

Suppose our 340-foot shot is hit down the line where the outfield fence is 345 feet away, so it's caught just in front of the wall. Assuming the initial velocity is still 160 ft/sec, the question is an what inclination angle – presumably somewhat greater than 30( -- we should hit it for a distance of 350 feet, and hence convert an out to a home run. The following function gives the distance hit, minus 350, as a function of the input angle alpha in degrees.

type fna

function distance = fna(alpha)

% distance = fna(alpha)

% input alpha in DEGREES

% outputs the (distance - 350) traveled if ball hit

% at angle alpha (with v0 = 160, c = 0.0025, k = 2)

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],160,alpha*pi/180,0.0025);

distance = ye(2,1)- 350;

So we seek a zero of this function. Note that

[fna(30), fna(40)]

ans =

-9.4491 2.8781

so there's a zero between alpha = 30( and alpha = 40(. We'll use MATLAB's fzero function to find it.

nflops = flops; tic;

alpha = fzero('fna', [30 40])

time = toc

nflops = flops - nflops

alpha =

34.3629

time =

14.6100

nflops =

332002

So an initial angle of about 34.36( will do it. This took about 14 seconds with a 266 MHz Pentium II machine cranking out

nflops/time

ans =

2.2724e+004

that is, about 23000 floating point operations per second, but the difference

fna(alpha)

ans =

-4.5475e-013

between 340.55 feet and our actual range is pretty close to a zero!

Fitting a Linear Velocity-Proportional Resistance Model

Suppose we ask what the proportionality constant c should be in order that with k = 1 (so air resistance R = cv) our standard initial conditions v0 = 160 ft/sec and ( = 30( would yield the same horizontal range 340.55 feet as we found with velocity-squared resistance R = 0.0025 v2. The following function of c measures the discrepancy between the actual range and the target value 340.55.

type fnc

function distance = fnc(c)

% distance = fnc(c)

% input the coefficient c in R = cv

% outputs the (distance - 340.55) traveled if ball hit

% at angle alpha (with v0 = 160, c = 0.0025, k = 1)

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],c,1);

distance = ye(2,1)- 340.55;

We note that

[fnc(0.2) fnc(0.4)]

ans =

63.2337 -64.5267

so there must be a zero near c = 0.3:

c = fzero('fnc', [0.2 0.3])

c =

0.2823

Thus linear air resistance R = 0.2823 v (with the default initial conditions) should yield the same range of 340.55 feet as does velocity-squared resistance R = 0.0025 v2. Let's check it out:

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],c,1);

te

te =

1.8917

4.1920

ye

ye =

203.0959 68.9658 81.2383 0

340.5500 0.0000 42.4406 -54.1436

Note that we have a trajectory very similar to the case of velocity-square resistance, except with a maximum height of 68.97 feet (rather than a maximum height of 65.64 feet).

Fitting a More General Air Resistance Model

Suppose a ball is hit with v0 = 160 ft/sec, alpha = 30( -- and careful radar measurement of the resulting trajectory shows a range of 340.55 feet but a maximum height of exactly 68.5 feet, between the max heights of 65.64 feet and 68.97 feet obtained with k = 1 and with k = 2, respectively. We then suspect air resistance like R = c vk with 1 < k < 2. We want to find both c

and k.

Given c and k, the following function calculates the sum of the squares of the differences between the range and max height and their target values.

type sqerr

function sse = sqerr(p)

% sse = sqerr(v)

% input the parameter vector p = [c k]

% outputs the sum of squares of "errors"

% resulting from the model R = c v^k

c = p(1);

k = p(2);

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],c,k);

sse = (ye(1,2)- 68.5)^2 + (ye(2,1)- 340.55)^2;

Obviously we want to find a zero of this function of two variables – which should be the same as a minimum value. But before we can use the MATLAB function fmins, we need a fairly good initial guess. So let's start with some trial-and-error, thinking that perhaps k ( 1.5.

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],0.1,1.5); ye

ye =

98.2423 38.5655 44.3163 0.0000

147.1931 0.0000 14.6754 -35.5634

Obviously k = 1.5 is too large – the ball goes nowhere near 340 feet, because there's too much air resistance. So let's decrease the value of k gradually until we're in the ball park (so to speak).

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],0.1,1.4); ye

ye =

131.0324 48.9106 55.8542 0.0000

203.9564 0.0000 21.9704 -42.0808

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],0.1,1.3); ye

ye =

168.9882 59.7838 69.5740 0.0000

275.4743 0.0000 33.3235 -49.2158

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],0.1,1.2); ye

ye =

208.7849 70.1372 84.4111 0.0000

357.4210 0.0000 49.1823 -56.4560

So k = 1.2 looks like a fairly good guess, but its still going a little too high and a little too far, soperhaps we should increase c a bit.

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],0.11,1.2); ye

ye =

201.0045 68.2109 81.3820 0.0000

340.4239 0.0000 45.3717 -54.9468

This is pretty close to the target of ymax = 68.5, xmax = 340.55. So now let's use fmins to solve the equation sqerr(c,k) = 0.

nflops = flops; tic;

v = fmins('sqerr', [0.11 1.2])

time = toc

nflops = flops - nflops

v =

0.1102 1.1995

time =

35.7600

nflops =

811766

Over 0.8 million floating point operations – for very little change in the input values. So we wonder whether the result is worth the effort?

sqerr(v)

ans =

0.0736

c = v(1); k = v(2);

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],c,k);

ye

ye =

201.0872 68.2321 81.4118 0.0000

340.5929 0.0000 45.4022 -54.9602

We're still off in the 2nd decimal place of each targeted value. The default MATLAB function fmins uses a relatively unsophisticated Nelder-Mead type simplex search method. We could get closer by decreasing the RelTol in bbode3, but instead let's try the much more sophisticated leastsq function in the MATLAB Optimization Toolbox. The following function defines the vector of individual errors that leastsq needs.

type errors

function errs = errors(v)

% errs = errors(v)

% input the vector v = [c k]

% outputs vector of two "errors"

% resulting from the model R = c v^k

c = v(1);

k = v(2);

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],c,k);

errs = [(ye(1,2)- 68.5); (ye(2,1)- 340.55)];

nflops = flops; time = clock;

v = leastsq('errors', [0.1 1.2])

time = etime(clock, time)

nflops = flops - nflops

v =

0.1570 1.1244

time =

75.3000

nflops =

1988455

Now almost 2 million floating point operations at

nflops/time

ans =

2.6407e+004

that is, over 26 thousand flops. To check the accuracy:

c = v(1)

k = v(2)

[T,Y,te,ye,ie] = ode45('bbode3',[],[],[],[],[],c,k);

ye

c =

0.1570

k =

1.1244

ye =

201.8082 68.5000 81.3349 0.0000

340.5500 0.0000 44.3205 -54.6750

So the resulting air resistance model

[pic].

looks pretty good – at least 4-decimal place accuracy in each target value. In fact, we have solid 5-place accuracy:

errors(v)

ans =

1.0e-006 *

-0.0662

0.1965

sqerr(v)

ans =

4.3001e-014

Not bad.

Copyright C. Henry Edwards 1999 Athens GA 30606

................
................

In order to avoid copyright disputes, this page is only a partial summary.

Google Online Preview   Download