Challenge¶

In this reading, your challenge is to use what you have learned to predict how the velocity of the blue car will evolve during the following experiment, after it reaches the flat portion of the table. You may collect data from this video for model evaluation only.

To aid you in constructing a model that you can use for this task, Data were collected for a similar experiment using a high-accuracy sensor, and are provided with this assignment. Columns in the data file are time and car velocity (m/s). You can load the data using the following command:

data = load('reading6data.txt');

You will need to examine how the derivative of the car's velocity (its acceleration) relates to the velocity itself during the experiment. You will be able to derive a differential equation for the system if you develop a model for this relationship.

In the last assignment, you used repeated extrapolations to get a solution to a differential equation. You can do the same here, but to avoid the tedious nature of repeating the extrapolation calculations manually over and over, it will help to develop a generalizable algorithm. This "repeated extrapolation" you performed in the last assignment is actually a form of approximate integration. Because the differential equation for some function $Y$ involves that function's derivative or derivatives, integration is a method of finding a solution to the differential equation.

Integrals: Definition¶

For our purposes in this course, we will define the integral of a function $y(t)$, using the notation $\int y(t) dt$, as follows:

The function $Y(t)$ is the integral of a function $y(t)$ if $\frac{d}{dt} Y(t) = y(t)$ for all values of $t$.

This means, in simple terms, that the derivative of $Y(t)$ is $y(t)$. Integrals are often visualized or introduced in calculus classes as "the area under a curve," which is certainly a useful way to conceptualize them. Consider, for example, the scenario in the graphic below, which explains the relationship between a constant and its integral.

When the functional form of $y$ is known, it is possible to compute an indefinite integral for it, yielding the function $Y$ in our definition. However, for an arbitrary dataset, whose functional form may not be known, often the best we can do is to approximate the integral of a function.

The process of "repeated extrapolation" that you used in the previous assignment is actually a form of approximate integration. How? Read the definition of approximate integration using "Euler's Method" below.

Approximate Integration: Euler's Method¶

In Euler's method of approximate integration, also called the "rectangular approximation" of an integral, the area under a curve, $y(t)$, which yields $\int y(t) dt = Y(t)$, can be approximated by repeatedly summing the area of rectangles projected either forwards or backwards in time from each discrete measurement of $y$. The process begins with some initial, known value of the integral $Y_0$, the value of $Y(t)$ is incremented by the area of the corresponding rectangle on the plot of $y(t)$. The process is shown visually below.

As you can see in the second plot, the idea is to use one, known value of our integral, $Y$, at the very start of the procedure. This is often called the initial value of $Y$. Then, the area under the derivative plot from $y(t_1)$ to $y(t_2)$ is added to our estimate of $Y$, and this value, when added to $Y_0$, becomes $Y(t_2)$. Formally, this can be written as:

$Y(t_2) = Y_0 + y(t_1)\cdot(t_2-t_1)$

Aside: Compare this to what you know about point-slope form for a line: $y = m(x-x_0) + y_0$, with $x_0$ and $y_0$ describing a known point on the $y$ (dependent variable) vs. $x$ (independent variable) plot. It's the same equation!

From $t_2$ to $t_3$, the procedure is repeated. In other words, the rectangle defined by the value of $y(t_2)$ is used to add "a little more" to our estimate to produce $Y(t_3)$. Formally,

$Y(t_3) = Y_0 + y(t_1)\cdot(t_2-t_1) + y(t_2)\cdot(t_3-t_2)$

Using a "recursive formulation," we could equivalently write:

$Y(t_3) = Y(t_2) + y(t_2)\cdot(t_3-t_2)$

And so on. In this description, we projected the value of $y$ forwards in the creation of our approximate rectangles. Therefore, this particular variant of Euler integration is called "Forward Euler" integration. "Backward Euler" integration is also possible. As with approximations of derivatives, the smaller the differences between times we use in the projection of the rectangles, the better our approximation of the derivative will be.

Because approximate integration is often performed on data arranged as arrays, we can use the counting index $k$ to write this as an equation that could be built into a "for loop" and computed efficiently in code. Considering arrays representing time $t$, derivative $y$, and integral $Y$, we could write an equation for any entry in the output array $Y(k)$ as:

$Y(k) = Y(k-1) + y(k-1)\cdot(t(k)-t(k-1))$

This should look similar to the equation for a line that you might use for linear extrapolation. Indeed, this equation is an equation for a line in point-slope form. We use each estimate for $Y(k-1)$ to produce an estimate for the next value of our integral $Y(k)$ using linear extrapolation. This is shown visually on the right panel of the plot above.

An example of how to implement approximate Euler integration in MATLAB/Octave is shown below, given a vector of known values for $y$. The approximate integral is compared with the "true" integral of the function to show its accuracy. Try changing the timestep $T$ in the code cell below to see how accuracy is affected by the size of the time step used for extrapolation/integration. A large value of $T$ is problematic for the same reason large extrapolations were problematic in Reading 4.

T = 0.01; %timestep
t = 0:T:1; %time vector
y = cos(10*t); %function we wish to integrate. This may come from data rather than a known function of t.

Y_estimate = zeros(size(t)); %initialize an array to hold our estimates for Y, the integral of y.
Y0 = 0; %some initial, known value for the integral at the beginning.
Y_estimate(1) = Y0; %place this value in the first entry of our array

Y_true = 0.1*sin(10*t); %for this example, we know that d/dt(sin(t)) = cos(t). 
%We will use this "exact" solution to the integral to allow us to see how accurate Euler integration is.

%now set up a for-loop starting at the second entry in Y_estimate, since we already know Y_estimate(1)=Y0
for k=2:length(t) %automatically set up the for-loop to be the right length
    %now implement Euler's formula. remember that here, y is dY/dt.
    Y_estimate(k) = Y_estimate(k-1) + y(k-1)*(t(k)-t(k-1));
end

%now we will compare our approximate integral with our "true" value to see if Euler integration works.

figure
plot(t,Y_true,'k.',t,Y_estimate,'r')
xlabel('Time (s)')
ylabel('Y(t)')
legend('true','Forward Euler Integration')

As you can see, Euler integration can be fairly accurate. Try different values of the time step $T$ in the code cell above to see how its value affects the accuracy of Euler integration.

Euler's method for differential equations¶

Introduction¶

To perform Euler integration to solve a differential equation, the process is exactly the same, except that we do not have a vector of known values for $y$! We will instead use the differential equation to get the "slope estimates" we need for numerical integration.

Imagine that you produced a model $\frac{dy}{dt} = -4y +2$, where your dependent variable was $y$ and your independent variables were the parameters "4" and "2," along with time $t$.

Feedback in differential equations¶

As we can see, $\frac{dy}{dt}$ depends on the value of $y$ itself, which makes Euler's method slightly more complex. How? If we were to visualize this relationship using a special flow chart called a "block diagram," we might draw: The block diagram is a visual representation of the fact that information about the value of $y$ is "fed back" into the computation of $\frac{dy}{dt}$. The round circle represents addition, and rectangles represent "operations" like multiplication and integration. This means that we can see that the equation: $$\frac{dy}{dt} = -4y + 2$$ Is represented in the diagram, with the result of the addition of $-4y$ and $2$, which is $\frac{dy}{dt}$, being sent into an integration function.

Euler's method for differential equations: implementation¶

The feedback in our differential equation means that in our "for loop" for integration, we have to use our most recent estimate for $y$ when computing our derivative $\frac{dy}{dt}$. This has to happen at every time step! See the example below, which uses Euler's method to solve our differential equation.

clear all
%create a time vector
t = 0:.001:1;
%create a vector of zeros to hold our y estimates
y = zeros(size(t));
%declare the first element in "y" as our best guess for the FIRST value of y (could be from a dataset, or knowledge of an experiment)
y(1) = 0.2;

%loop through time values to perform integration. start at the second value of time, since we already have a guess for y at the first value!
for k=2:length(t)
    %first use the differential equation to compute dy/dt (slope) based on our best guess for y at the previous time
    dydt = -4*y(k-1) + 2; %2 is a constant, but it doesn't have to be. I could have a vector of known values for 2 as well.
    %now extrapolate (integrate) to our current time to get a new estimate for y!
    y(k) = dydt*(t(k)-t(k-1)) + y(k-1);
end

figure
plot(t,y,'k')
xlabel('Time (s)')
ylabel('y')

As you can see, Euler's method performs the same "repeated extrapolation" you did by hand in Reading 5 without knowing its formal name, this time with a convenient for-loop!

In your challenge for this week, you may find Euler's method, which is a quick way to perform the 'repeated extrapolation' you used in Reading 05, helpful.

Note that in order to use Euler's method, you do have to have a differential equation!. In Reading 05, you made a crude guess for this differential equation based on data, but it is possible to be a little bit more precise when deriving a data-driven, empirical differential equation. We can be more precise by fitting a relationship between $\frac{dy}{dt}$ as the independent variable and $y$ as the dependent variable.

Empirical Models: fitting polynomials in MATLAB/Octave¶

MATLAB/Octave can be used to estimate a "best fit" polynomial function to match a dataset. The term "best fit," in this context, means that the sum of the squared error between the model and the data is minimized. The command used to produce the set of coefficients is "polyfit." An example of how to use the "polyfit" command to produce the best-fit line for a dataset is shown below.

% first, create some fake "data" we will fit a line to, in order to see how "polyfit" works.
xfakedata = 0:.1:10;
yfakedata = -3*xfakedata + 4 + .5*randn(size(xfakedata)); %the randn() on the end produces some random "noise" on our fake data.
%this simulates the effect of having error in our measurements. If polyfit does its job, we should get the line y = -3x+4 back!

%now we will actually use polyfit. The array "coefficients" will contain [m b] in the equation y = m*x+b.
%the second argument to the polyfit() function is the ORDER of the polynomial. A line is a polynomial of order 1. 
%for reference, a quadratic is a polynomial of order 2.
coefficients = polyfit(xfakedata,yfakedata,1); 
%was m=-3 and b=4?
m_estimate = coefficients(1)
b_estimate = coefficients(2)

%now, these coefficients fully define a "model" of our relationship y = m*x+b. Let's produce predictions based on this model
ymodel = m_estimate*xfakedata + b_estimate;

%now plot the data and the model together to evaluate.
figure
plot(xfakedata,yfakedata,'ks',xfakedata,ymodel,'r')
xlabel('x')
ylabel('y')
legend('data','model')

m_estimate = -2.9648
b_estimate =  3.9016

Assignment¶

Using your now-refined knowledge of repeated extrapolation, or the "Forward Euler" approximation for integration, along with your newfound knowledge of how to estimate a functional relationship between two variables using the "polyfit" command, you will:

Use the data provided in reading6data.txt to estimate a functional relationship between the car's velocity and its acceleration after it reaches flat ground
Use this functional relationship, along with the first measured speed of the car in the video provided, to predict the car's velocity during the experiment in the video in 0.01 second increments. Note that you will have to create your own time vector, since the video's frame rate is not that high.
Use this prediction of the car's velocity during the experiment to estimate its position during the experiment. You can obtain position estimates by integrating estimates your velocity model using Euler's method.
Using measurements of car position from the video above, compare your model with data from the experiment to evaluate your model for the car's position.

To make data collection for the model evaluation step easier, a slowed-down and timestamped version of the original video is provided below. Rapid use of the spacebar should allow you to record a position of the car for every unique time in the video.

Deliverables: System Scoping¶

Provide a diagram showing your system model's inputs and outputs, along with a terse but complete description of the contents of the model.

YOUR ANSWER HERE

Deliverables: Model Construction¶

Provide any hand calculations you use in the construction of your model in the Markdown cell provided, and any code you use to fit a functional relationship between the car's velocity and its acceleration in the code cell below. You may also "run your model" in this code cell if you wish to provide estimates of the car's velocity during the test. The dataset in reading6data.txt has been provided for you to use in your model construction steps. It provides two columns. The first column is time in seconds, and the second column is velocity in meters/second.