% example mixing problem % by Gavin LaRose, 6/28/2012 % this demonstration models the number of black balls in a box as % a mixing problem. we use a continuous model and plot the % solution to the differential equation. % % number of black balls added each time unit blackInPerFrame = 3; % the size of the hole through which balls leave the box, in % ball units holeSize = 4; % how long to calculate for simulationLength = 200; % initial number of black balls in the box numBlack = 0; % "volume" of the box, in number of balls; should be divisible % by the boxWidth, below. good values are 50, 100 and 200 volume = 100;
The code block above just sets up the variables describing the problem.
% (numerical) solution to the modeling differential equation; this
% gives the number of black balls in the box
[t,theorySol] = ode45(@(t,y) blackInPerFrame - holeSize*y/volume, ...
1:simulationLength, numBlack);
The arguments for ode45 are, in order, a function handle
giving the value of \(y'\) in the equation we are solving, the time range
over which to solve, and the initial condition. The @ sign
indicates that we have a function handle, and in this case we declare it
as an anonymous function of two variables (t,y); we have \(y'
= \mbox{input} - h\,y/\mbox{volume}\), so the function returns
blackInFrame - holeSize*y/volume. We start the time at
\(t=1\) rather than \(t=0\) so that it matches the simulation that we do
in Example_Mixing_Discrete.m and
Example_Mixing_Problem.m. Note that the construction
1:simulationLength creates a vector of values, \(\lt 1, 2,
3,\ldots\ \mbox{simulationLength}\gt\).
Obviously, we could instead give the exact solution to this differential equation, but the (approximate) numerical solution will be fine for our purposes.
% we can then plot this result
figure();
plot(t, 100*theorySol/volume, 'b-','LineWidth',2);
axis( [0 simulationLength 0 volume] );
title_handle=title('Balls in a Box: ODE Model');
set( title_handle, 'FontSize', 18 );
legend_handle=legend( '% of black balls in box' );
set( legend_handle, 'FontSize', 14, 'Location', 'SouthEast' );
In the plot, we give vectors for the \(x\)- and \(y\)-coordinates of the
points to plot, then specify the plot type (b- is blue, with
solid line), and increase the line width to make it more visible. Note
that in the calculation of the percent of balls in the box we are taking
advantage of Matlab's vector operators: 100*theorySol/volume
takes \(100/\mbox{volume}\) times each element of the vector
theorySol. The axis command sets the length of
the \(x\)- and \(y\)-axes, in order. We save the handles of the title and
legend objects so that we are able to change their font size and set other
characteristics for them.