- 42 Views
- Uploaded on
- Presentation posted in: General

Bruce Mayer, PE Registered Electrical & Mechanical Engineer [email protected]

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -

Engr/Math/Physics 25

Chp6 MATLABFcn Discovery

Bruce Mayer, PE

Registered Electrical & Mechanical [email protected]

- Create “Linear-Transform” Math Models for measured Physical Data
- Linear Function → No Xform
- Power Function → Log-Log Xform
- Exponential Function → SemiLogXform

- Build Math Models for Physical Data using “nth” Degree Polynomials
- Use MATLAB’s “Basic Fitting” Utility to find Math models for Plotted Data

- Use Regression Analysis as quantified by the “Least Squares” Method
- Calculate
- Sum-of-Squared Errors (SSE or J)
- The Squared Errors are Called “Residuals”

- “Best Fit” Coefficients
- Sum-of-Squares About the Mean (SSM or S)
- Coefficient of Determination (r2)

- Sum-of-Squared Errors (SSE or J)
- Scale Data if Needed
- Creates more meaningful spacing

- Calculate

- Build Math Models for Physical Data using “nth” Degree Polynomials
- Use MATLAB’s “Basic Fitting” Utility to find Math models for Plotted Data
- Use MATLAB to Produce 3-Dimensional Plots, including
- Surface Plots
- Contour Plots

Physical Processes for some Response (OutPut), y, as Resulting from some Excitation (InPut), x, can many times be approximated by 3 Functions

The LINEAR function: y = mx + b

Produces a straight line with SLOPE a of m and an INTERCEPT of b when plotted on rectilinear axes

The POWER function: y = bxm

gives a straight line when plotted on log-log axes

The exponential function y = b·10mx or y = b·emx

Yields a straight line when plotted on a semilog plot with logarithmic y-axis

For a Linear Function We can Easily Find the Slope, m, and y-Intercept, b

- Transform the Power Function to Line-Like Form

- Take ln of both sides

Thus the Power Function Takes the Form:

- Yields a Staight Line
- So if we suspect a PwrFcn, Plot the Data in Log-Log Form

- Example:

- If The Data Follows the Power Function Then a Plot of

Rectlinear Plot

Log-Log Plot

>> x = linspace(7, 91, 500);

>> y = 13*x.^1.73;

>> loglog(x,y), xlabel('x'), ylabel('y'), title('y = 13x^1^.^7^3'), grid

Recall the General form of the Exponential Fcn

- In This Case Let

- Then the Xformed Exponential Fcn

- Again taking the ln

- A SemiLog (log-y) plot Should Show as a Straight Line

Plot

- SemiLog Plot

- Rectlinear Plot

>> t = linspace(0, 10, 500);

>> v = 115*exp(-0.61*t);semilogy(t,v), xlabel('t'), ylabel('v'), title('v = 115e^-^0^.^6^1^t'), grid

Examine the data near the origin.

The linear function can pass through the origin only if b = 0

The exponential function (y = bemx) can never pass through the origin

as et > 0 (positive) for ALL t; e.g., e−2.7 = 0.0672

unless of course b = 0, which is a trivial case: y = 0·emx

The power function (y = bxm) can pass through the origin (e.g.; y = 7x3) but only if m > 0 (positive) as

As y = bx-m = b/xm→ Hyperbolic for negative m

- In most applications x is NONnegative

Plot the data using rectilinear scales.

If it forms a straight line, then it can be represented by the linear functionand you are finished.

Otherwise, if you have data at x = 0, then

If y(0) = 0, then try the power function.

If y(0) 0, then try the exponential function

If data is not given for x = 0, then proceed to step 3.

If you suspect a power function, then plot the data using log-log scales.

Only a power function will form a straight line on a log-log plot.

If you suspect an exponential function, then plot the data using the SemiLogy scale.

Only an exponential function will form a straight line on a SemiLog plot.

- Note Change in power Function x-axis scales
- In this case x MUST be POSITIVE

In function discovery applications, use the log-log and semilog plots only to identify the function type, but not to find the coefficients b and m.

The reason is that it is difficult to interpolate on log scales

To Determine Quantities for m & b Perform the Appropriate Linearization Transform to plot one of

ln(y) vs. ln(x) → Power Fcn

ln(y) vs x → Exponential Fcn

The Command → p = polyfit(x,y,n)

This Function Fits a polynomial of degree n to data described by the vectors x and y, where x is the independent variable.

polyfit Returns a row vector p of length n + 1 that contains the polynomial coefficients in order of descending powers

Note That a FIRST Degree Polynomial Describes the Eqn of a LINE

If w =p1z +p2 then polyfit on data Vectors W & Z returns: p = ployfit(Z,W,1) = [p1, p2] → [m, b]

polyfit of degree-1 (n = 1) returns the parameters of a Line

p1 → m (slope)

p2 → b (intercept)

Thus polyfit can provide m & b for any of the previously noted functions AFTER the appropriate Linearization Transform

We need to find an Eqn for the Vapor pressure of Ethanol, C2H5OH, as a fcn of Temperature

Find Pv vs T Data by Consulting

P. E. Liley and W. R. Gambill, Chemical Engineers’ HandBk, New York, McGraw-Hill Inc., 1973, p. 3-34 & 3-54

Since this a Vapor Pressure we Suspect an Antoine or Clapeyron Relation Pv ~ em/T

As a Starting Point Make a Rectilinear Plot

>> Pv = [1, 5, 10, 20, 40, 60, 100, 200, 400, 760];

>> T = [241.7, 261, 270.7, 281, 292, 299, 307.9, 321.4, 336.5, 351.4];

>> plot(T,Pv,'x', T,Pv,':'), xlabel('T (K)'), ylabel('Pv (Torr)'),...

title('Ethanol Vapor Pressure'), grid

Log Log Plot

- logY vs linX

- Looks like an Exponential in 1/T
- e.g., the Clapeyron eqn

The Plots Looks Pretty Well Exponential

For a 1st-Cut Assume the Clapeyron form

- Now Xform

- Thus Plot ln(Pv) vs 1/T

The Command Session

>> Ye = log(Pv);

>> x = 1./T

>> plot(x,Ye,'d', x,Ye,':'), xlabel('1/T (1/K)'), ylabel('ln(Pv) (ln(Torr))'),...

title('Ethanol Vapor Pressure - Clapeyron Plot'), grid

- Nicely Linear → Clapeyron is OK

Apply PolyFit to Find m and B

- To Increase the Sig Figs displayed for B these plots are typically plotted with x1 = 1000/T

>> p = polyfit(x,Ye,1)

p =

1.0e+003 *

-5.1304 0.0213

>> x1 = 1000./T

>> p1 = polyfit(x1,Ye,1)

p1 =

-5.1304 21.2512

- or

PowerandExps

Engr/Math/Physics 25

Appendix

Time For

Live Demo

Bruce Mayer, PE

Licensed Electrical & Mechanical [email protected]

>> x =[7:91];

>> y = 13*x.^1.73;

>> Xp = log(x); >> Yp = log(y);

>> plot(x,y), xlabel('x'), ylabel('y'),...

grid, title('Power Function')

>> plot(Xp,Yp), xlabel('ln(x)'), ylabel('ln(y)'),...

grid, title('Power Function')

>> t = [0:0.05:10];

>> v = 120*exp(-0.61*t);

>> plot(t,v), xlabel('t'), ylabel('v'),...

grid, title('Exponential Function')

>> Ve = log(v);

>> plot(t,Ve), xlabel('t'), ylabel('ln(v)'),...

grid, title('Exponential Function')

The Area of RIGHT Triangle

- The Area of an ARBITRARY Triangle

- By Pythagoras for Rt-Triangle

- Equating the A=½·Base·Hgt noting that

- Solving for h

Q.E.D.

>> T = [69.4, 69.7, 71.6, 75.2, 76.3, 78.6, 80.6, 80.6, 82, 82.6, 83.3, 83.5, 84.3, 88.6, 93.3];

>> CPH = [15.4, 14.7, 16, 15.5, 14.1, 15, 17.1, 16, 17.1, 17.2, 16.2, 17, 18.4, 20, 19.8];

>> Tmax = max(T);

>> Tmin = min(T);

>> CPHmax = max(CPH);

>> CPHmin = min(CPH);

>> Rtemp = Tmax - Tmin;

>> Rcroak = CPHmax - CPHmin;

>> DelT = T - Tmin;

>> DelCPH = CPH - CPHmin;

>> Theta = DelT/Rtemp;DelCPH = CPH - CPHmin;

>> Omega = DelCPH/Rcroak;

>> plot(T, CPH,), grid

>> plot(Theta,Omega), grid

FIRST → Plot the DATA

Goodness of Fit; smaller is Better

Expand Dialog Box

Result

Chk by polyfit

>> p = polyfit(Theta,Omega,1)

p =

0.8737 0.0429

% "Discoverable" Functions Displayed

% Bruce Mayer, PE • ENGR25 • 15Jul09

%

x = linspace(-5, 5);

ye = exp(x);

ypp = x.^2;

ypm = x.^(-2);

% plot all 3 on a single graphe

plot(x,ye, x,ypp, x,ypm),grid,legend('ye', 'ypp', 'ypm')

disp('Showing MultiGraph Plot - Hit ANY KEY to continue')

pause

%

% PLot Side-by-Side

subplot(1,3,1)

plot(x,ye), grid

subplot(1,3,2)

plot(x,ypp), grid

subplot(1,3,3)

plot(x,ypm), grid