14 Numerical Computing

Since this book (and the course it is based on) focuses on programming, it has not emphasized MATLAB’s many built-in capabilities for performing mathematics. In reality, it is almost always smarter to use “canned” routines for performing complex computations, since many hours of development have gone into ensuring their accuracy and efficiency. While writing a simulation from scratch using Euler’s method is a valuable educational exercise, it is not the best use of a practicing engineer’s time, when MATLAB has functions to perform the calculation with little or no programming.

This chapter covers some of MATLAB’s built-in functions for several kinds of fundamental engineering analysis:

• curve fitting
• interpolation
• solving a system of linear equations
• solving an ordinary differential equation
• optimization and parameter estimation

14.1 Polynomial Curve Fitting

Extracting information from experimental data is a fundamental task in engineering analysis. A basic technique for accomplishing that is to fit a curve describing a theoretical relationship to the data. The standard approach is least-squares fitting, in which the model parameters (e.g. the coefficients of a polynomial) are chosen to minimize the sum of the squares of the differences between the theoretical curve and the experimental points. MATLAB has a built-in function, polyfit(), to perform least-squares fitting with polynomial models. The input arguments are the x and y data arrays and the order, n, of the polynomial. The output argument is an array containing the coefficients of the polyomial.

As a first example, consider a linear fit to a set of deflection vs. load data for a cantilever beam and extraction from the elastic modulus from the slope of the line. The theoretical relationship between deflection, [latex]\Delta y[/latex], and the applied load, P, for an end-loaded cantilever beam having an elastic modulus E is:

[latex]\Delta y = \frac{PL^3}{3EI}[/latex]

where L is the free length of the beam, and I is the area modulus of elasticity, which depends on the cross-sectional area and shape of the beam. For a solid rectangular beam of width w and height h

[latex]I = \frac{1}{12}wh^3[/latex]

A common first-year engineering experiment is to measure the deflection of a beam for several applied loads. A plot of [latex]\Delta y[/latex] vs. P should be a straight line with slope

[latex]m = \frac{L^3}{3EI}[/latex]

Suppose the experimental data are as given in arrays delta_y and P. The coefficients of the least-squares line are calculated by the polyfit() function as follows:

coeffs = polyfit(P, delta_y, 1);

The third argument, 1, indicates that it is a first-order (linear) polynomial. The function returns a vector of coefficients in descending order; in this case coeffs(1) is the slope, and coeffs(2) is the y-intercept.

Once the coefficients of the fitted linear are determined, a plot of that line can be created with the help of the polyval() function. This function takes as inputs the vector containing the coefficients and an array of x-values at which the corresponding y-values are to be calculated. It returns the vector of y-values, which can be used to plot the curve. The example below shows the complete process of fitting a straight line, generating a plot of the data overlaid with the best-fit line, and calculating the elastic modulus based on the slope.

14.2 Interpolation

A closely related topic to curve fitting is interpolation, which means estimating the value of a quantity between known values. For example, a transducer – let’s say a pressure sensor – may have been calibrated by making measurements at several known pressure values. Suppose that the calibration data is as follows:

After calibration, an unknown pressure is measured, and the resulting voltage is 56.5 mV. What is the pressure? It must be between 1.0 and 1.1 atm, since the measured voltage is between the voltage measurements at those two calibration points, but what precisely? That depends on the specific interpolation method used; the simplest method is linear interpolation, in which the dependent variable is assumed to vary linearly with the independent variable between data points. See the linked article for details of the calculation.

MATLAB has a built-in function, interp1(), to perform one-dimensional interpolation. It accepts as arguments:

• a vector containing the known x-values
• a vector containing the known y-values
• a vector containing the x-values at which the y-values are to be estimated

For the example of the pressure sensor described above, the function call could be:

>> P = interp1([52.0, 68.3, 79.5], [1.0, 1.1, 1.2], 56.5)

P =

1.0276

>>

The interp1() function can perform other interpolation methods (see the documentation for details). One commonly used method is cubic spline interpolation, which draws a smooth curve through the data points, rather than connecting them by straight lines. To choose spline interpolation instead of linear, an additional argument is passed to interp1() as follows:

>> P = interp1([52.0, 68.3, 79.5], [1.0, 1.1, 1.2], 56.5, 'spline')

P =

1.0222

>>

The example below performs both linear and spline interpolation on a set of data and plots them to illustrate the difference.

Interpolation vs. Curve Fitting

14.3 Systems of Linear Equations

14.4 Ordinary Differential Equations

14.5 Optimization

Find an error? Have a suggestion for improvement? Please submit this survey.