We discuss Modeling and Parameter Estimation

One of the most common ways that the mathematical modeling structure can be used to analyze physical problems is the idea of parameter estimation. The situation is that we have physical principles that give rise to a differential equation that defines how a physical system should behave, but there are one or more constants in the problem that we do not know. Two simpler examples of this are Newton’s Law of Cooling

which models the temperature of an object in an environment of temperature over time, and velocity affected by drag modeling the velocity of a falling object where the drag force is proportional to the square of the velocity. In both of these cases, the models are well established, but for a given object, we likely do not know the or values in the problem. These are these parameters of the problem, and would be determined by the shape and structure of the objects, the material that it is made of, and many other factors, so it could be hard to figure out what they are in advance. How can we find these values? We can use data from the actual physical problem to try to estimate these parameters.

The easier version of this is to use a single value at a later time to calculate the constant.

Solution: Based on Newton’s Law of Cooling, we know that the temperature satisfies the differential equation

with initial condition , but we do not know the value of . In order to work this out, we should solve the differential equation with unknown constant , then figure out which value of gives us the appropriate temperature after 10 minutes. This is a first order linear equation, which can be rewritten as

The integrating factor we need is , which turns the equation into

Integrating both sides and solving for gives

To satisfy the initial condition, we need that , or . Thus, our solution, still with an unknown constant , is

To determine the value of , we need to utilize the other given piece of information: that . Plugging this in gives that

which we can solve for using logarithms. This will give that

Finally, we can plug that constant into our equation to get the solution for the temperature at any time value,

__

This method works great if we have the exact measurement from the object at one point in time. However, if the measurements at multiple points in time are known, and if the data is not likely to be exact, then a different method is more applicable. The idea is that we want to minimize the “error” between our predicted result and the physical data that we gather. The method used to minimize the error is the “Least Squared Error” method.

Assume that we want to do this for the drag coefficient problem,

where we do not know, and want to estimate, the value of . For this method, the data that we gather is a set of velocity values that are obtained at times . For any given value of , we can solve, either numerically or analytically, the solution to the given differential equation with that value of . From this solution, we can compute , the value of this solution at each of the times that we gathered data originally. Now, we want to compute the error that we made in choosing this parameter . This is computed by which is the sum of the squares of the differences between the gathered data and the predicted solution. In order to find the best possible value of , we want to minimze this error by choosing different values of and whatever value of gives us this minimum is the optimal choice for that parameter.

The function that we want to minimize here is usually a very complicated function, and we may not even be able to solve the differential equation analytically for any . Thus, computers are used most often here to solve these types of problems.

Solution: To solve this problem, we will use the least squared error method implemented in MATLAB. The code we need for this is the following, which makes use of the Optimization Toolbox.

1 global tVals
2 global vVals
3
4 tVals = [0,0.1,0.3,0.5,0.8,0.9,1.0,1.5];
5 vVals = [0,0.9797,2.8625,4.4750,6.3828,6.8360,7.0334,8.1612];
6
7 [aVal,errVal] = fminbnd(@(a) EstSqError(a),0,4)

This bit of code inputs the necessary values and uses the fminbnd function to find the minimum of the error function on a defined interval. These problems need to be done on a bounded interval, but in most physical situations there is some reasonable window for where the parameter could be. The rest of the code is the definition of the EstSqError function.

1 function err = EstSqError(al)
2
3 global tVals
4 global vVals
5
6 fun = @(t,v) 9.8 -al.*v.^2;
7 sol = ode45(fun,[0,3],0);
8 vTest = deval(sol,tVals);
9
10 err = sum((vVals -vTest).^2)
11 end

The main point of this code is that it takes in a value of , over which we are trying to minimize, it numerically solves the differential equation with that value of over a desired range of values, and then compares the inputted vVals with the generated vTest array, computing the sum of squared errors, and returning the error value.

Running this code results in an value of , with an error of . That means that, based on this data, the best approximation to is . ___

Note that in the above example, the total error was not zero, and doesn’t actually match the coefficient used to generate the data, which was . This is because noise was added to the data before trying to compute the drag coefficient. In a real world problem, noise would not be added, but a similar effect would arise from slightly inaccurate measurements or round-off errors in the data. While we may not have found the constant exactly, we got really close to it, and could use this as a starting point for further experiments and data validation.