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.
To fulfill the demand for quickly locating and searching documents.
It is intelligent file search solution for home and business.
Related download
Related searches
- minecraft apk v1.8.0
- lego education wedo software v1 2
- wedo software v1 2
- what s the best notebook computer
- windows system32 windowspowershell v1 0 powershell exe
- 5 subject notebook with pockets
- the truth project lesson 5 part b
- powershell v1 0
- system 32 powershell v1 0
- windows system32 windowspowershell v1 0 powershell
- evaluate the expression for b 5 2 1b
- icd 705 v1 4