Engineering Analysis EGN 3420
University of Central Florida
Popular in Course
Popular in General Engineering
This 146 page Class Notes was uploaded by Lora Metz on Thursday October 22, 2015. The Class Notes belongs to EGN 3420 at University of Central Florida taught by Staff in Fall. Since its upload, it has received 33 views. For similar materials see /class/227560/egn-3420-university-of-central-florida in General Engineering at University of Central Florida.
Reviews for Engineering Analysis
Report this Material
What is Karma?
Karma is the currency of StudySoup.
Date Created: 10/22/15
Section 4 Lagrange Interpolating Polynomials In the previous sections we encountered two different ways of representing the unique nth order or lower polynomial required to pass through a given set of n1 points Yet another way of writing the polynomial constrained in the same fashion is presented here It is referred to as Lagrange s form of the interpolating polynomial Once again we assume the existence of a set of data points x y 1 0 l 71 obtained from a function x so that y xi 1 0 l n A suitable function for interpolation 1x is expressible as n 1x ELImm 41 i0 L0x39fxoL1x39fx1 quotquot39Lnx39fxn 41a The functions Lix 1 0 l n are chosen to satisfy 0 x x0x1xi1xi1xn L x 42 1 x x 139 Before we actually define the Lix functions let39s be certain we understand the implications of Equations 41 and 42 The best way to accomplish this is simply to choose a value for quotnquot and write out the resulting equations Suppose we have the four data points xi xi 1 0 l 2 3 From Equations 41 with n 3 the interpolating function 1x becomes 3 1x lex r mi 43 i0 LOOC fx0L1x39fx1L2x39fx2L3x39fx3 43a and it remains to be shown that 1x is identical to x when x is any one of the four data points Evaluating 1x at x0 x1 x2 and x3 950 Lox039fxoL1x039fx1L2x039fx2L3xo39fx3 44 1x1 L0x139fxoL1x139fx1L2x139fx2L3x139fx3 443 1x2 L0x239fx0L1x239fx1L2x239fx2L3x239fx3 44b x3 L0x339fx0L1x339fx1L2x339fx2L3x339fx3 According to Equation 42 L0x0 l and L1x0 Lzx0 L3x0 0 Equation 44 is simplified as shown below x0 L0x039fx0L1xo39fx1L2xo39fx2L3xo39fx3 139fxo039fx1039fx2039fx3 fx0 By the same reasoning 1x1 xl 1xz xz and 1x3 xg which means the interpolating function Ix passes through the given set of data points The analytical form of Ix depends on the functions Lix 139 01 n that satisfy Equation 42 They are called Lagrange coefficient polynomials and are defined as follows x X39 Lix I I x i0l2n 45 j0li liln Z J The symbol H in Equation 45 is a symbolic notation for multiplication in the same way the symbol 2 denotes summation of its arguments To better understand Equation 45 Lx are expressed in a more explicit form Given below are expressions for L1x through L3x when there are four data points 71 3 w l 0 L000 x0 x1x0 x2x0 x3 46 11 Lx m 46a 1 x1 x0x1 x2x1 x3 I x XOXX xiXx x3 l 2 L200 962 xox2 x1x2 x3 4396b x XOXX xiXx xz l 3 L300 x3 x0x3 x1x3 x2 43960 Notice that each Lagrange coefficient polynomial in Equation 46 is a third order polynomial as a result of the x3 term in the numerator For the general case when there are n1 data points the Lagrange coefficient polynomials LX in Equations 41 are nth order polynomials and therefore so is the interpolating function 1x Henceforth we shall represent the interpolating polynomial Ix in Equation 41 by fnx and refer to it as the Lagrange interpolating polynomial The advantage of Lagrange interpolation in comparison with the standard polynomial form Equation 21 or the Newton divideddifference representation Equation 31 is its simplicity That is the Lagrange interpolating polynomial can be determined without need of solving a system of simultaneous equations or performing repetitive calculations as in the case of Newton interpolating polynomials where a table of divided differences is required Owing to the manner in which the Lagrange coefficient polynomials are defined in Equation 45 the Lagrange interpolating polynomial is written by inspection of the data points using Equation 41 We demonstrate the procedure in the following example Example 41 The measured voltage as a function of time across the terminals of a 5 ohm load with three different types of batteries each with a nominal rating of 15 volts is tabulated below t hours v volts Rechargeable Rechargeable Rechargeable NiCd Alkaline Alkaline 1St Charge 1St Charge 2quotd Charge 0 140 140 135 1 1 3 0 1 17 1 15 2 100 1 10 105 3 040 105 100 4 005 090 0 40 Table 41 Voltage and Elapsed Time for Rechargeable Batteries The Lagrange interpolating polynomial for each type of battery is obtained directly from the table For the NiCd battery it is tt1tt2tt3tt4 V0 tt0tt2tt3tt4 V0 to l 11 0 tzxto 1 3 to 1 4 0 t1 l 01 1 l 21 1 l 31 1 1 4 1 f4t tt0tt1tt3tt4 V0 tt0tt1tt2tt4 Vt3 t2 t0 t2 t1t2 t3t2 t4 t3 t0t3 t1t3 t2 t3 t4 tt0tt1tt2tt3 39VUA t4 tot4 t1t4 t2t4 t3 f I t lt 2t 3t 4 40 tt 2t 3t 4 130 0 10 20 30 4 39 1 01 21 31 4 39 tt 1t 3t 4 100 m 1t 2t 4 0 40 2 02 12 32 4 3 03 13 23 4 39 tt lt 2t 3 005 48 4 04 14 24 3 39 39 f4t t lt 2t 3t 4 15 311 2t 3t 4 12 0tt lt 3t 4 tt 1r 2t 4tt 1r 2t 3 49 The same procedure is used to obtain the fourth order Lagrange interpolating polynomials for the two alkaline batteries Figure 41 shows the data points and the fourth order interpolating polynomial for each battery It is generally advisable to begin with lower order polynomials for interpolation and then increase the order in a systematic fashion until the results are acceptable In other words a reduced subset of data points is chosen based on the expected range of the interpolant Suppose we have n1 data points measured from a function x where the range of x is from xmin to xmax Assuming interpolation over the entire interval is required with a single polynomial we might choose a third or possibly fourth order polynomial assuming ngt4 based on a judicious choice of four or five points that include both extremes Visual inspection of the polynomial and its relationship to the entire set of data points will determine if a higher order polynomial is required and which additional data point should be included When an additional data point is required the higher order polynomial will generally exhibit more curvature than the previous one throughout the interpolation interval For example Figure 42 shows polynomials of order n 1 2 3 and 4 required to pass through data sets of two three four and five points respectively A1ka1tne 1 charge A1ka1tne 2 d charge v vans n5 1 15 25 3 35 2 t huuts thure 4 1 Fourth order 1nterpo1aung Po1ynomta1s ForBatteryData tn Table 4 1 It s e1ear that the quadratte po1ynomta1 f2x exhtbtts more eurvature than the 1tnear funeuon xwh1h of eourse has none Ihts ts posstb1e due to the presenee of the x2 term tn stmt1ar1ysr curves more than2r due to the addtuona1 potnt x and the fouth order term r tn Q aeeounts for the added curvature eompared f3x overthe tnterva1 04 m to wtth an nth tnterpo1attng po1ynomta1 tn standard form Equauon 21 or the Newton ddvxdedrddfference form Equatton 3 1 there ts only one htgh order term t e a stng1e term wtth x Ihts ts tn eontrast to the Lagrange form of the tnterpo1attng po1ynomta1 Equatton 4 1 tn whteh eaeh term f e v expre st t th mta1 If 4 5 tnereases tn order from n to n1 as we11 Consequenuy the enttre Lagrange tnterpo1aung po1ynomta1 must be reeomputed Vudrdf vw mm w hth order term an extra data potnt requtres reea1eu1atton of a1 the coef ments a t r 1 n tn addtuon to the new eoetrtetent a The system of equauons to be solved was eonstdered tn Seetton 32 and enumerated tn matntr form tn Equauon 2 3 The Vandermonde mamx and the eo1umn veetors of coef ments and funeuon values tn Equauon 2 3 are n1gtlt n1ngtlt 1 anan 1 respeeuve1y The Newton ddvxdedrddfference rnterpolatrng polynomra1 of the three forms eonsrdered mrnrmrzes the eomputatrona1 etrort neeessary to aeeommodate an addmonal data pornt As we pornted out m Seeuon 3 3 the Newton ddvxdedrddfference rnterpolatrng polynomra1 x m Equatron 3 1 rs formulated m sueh away that the eoemerents 5 t V 39339 a u t ddfferences gt2 Frgure 4 2 Inereasrng curvature ofaPolynomxal as aFuneuon ofxts order th any rnterpolatron method based on the use of apprommauon funeuon there rs an error term present whreh aeeounts for the ddfference between the true funetron value usually unknown and the rnterpolated value Usmg an nth order Lagrange rnterpolaung polynomra1 x to estrmate values of x the error term Rm satrs es f0 A RX 410 whreh rs rdenuea1 to Equatron 3 31 m Seeuon 3 whenquotr3 was an nth order Newton r eneernterpolatrng poly omr nee mt ag n 1 d t pornts mu e an ysrs presented m Seeuon rs appheable to 3 agrng a1s hus an esumate of the e r 1w rs su1 the ddfference m x rx where an addmonal data pornt rs requrred to eva1uate x The real advantage of Newton divideddifference interpolating polynomials is the reduced computational effort compared with the Lagrange form to obtain 1x when x is already computed Example 42 In Example 41 estimate the error incurred by the use of a fourth order Lagrange interpolating polynomial for approximating the voltage of a NiCd battery after 15 hours Assume an additional data point is available namely t 25 hours v 07 volts From Equation 49 f4 15 15 115 215 315 4 1515 215 3X15 4 1505 115 315 4 1505 115 215 4 1515 115 215 3 f415 11965 Supplementing Table 41 with the new data point and determining the fth order Lagrange interpolating polynomial yields 1 1t 2t 3t 4t 25 W tt 2t 3t 4t 25 130 4 1 1 fs I 0 10 20 30 40 25 391 l 1 01 21 31 41 25 tt 1t 3t 4t 2 5 W I tt 1t 2r 4t 2 5 0 40 2 02 12 32 42 25 39 39 3 03 13 23 43 25 39 tt 1t 2r 3t 25 tt 1t 2r 3t 4 005 70 4 04 14 24 34 25 25 025 125 225 325 4 R415 is the difference between the true unknown function value f15 and the fourth order interpolating polynomial estimate f415 ie R 15 f15 f4 05 We can estimate R415 using R405 z 1215f415 rs 12148 11965 rs 00184 1 Exercises A simplified mathematical model of tra ic ow is based on the assumption that at any point along a road the velocity of a car depends only on the density of cars in the immediate vicinity At low densities drivers are free to drive at their desired speed called the mean free speed At high densities vehicle speeds are lower ultimately approaching zero at the maximum density corresponding to bumpertobumper traffic This model is applicable on single lane roads or tunnels where car passing is not allowed Measurements of average velocity in mph and traffic density vehiclesmile were taken in a tunnel to determine its velocitydensity pro le Density p Velocity u gvehicles mile 2 1 mph 2 205 398 312 309 504 227 864 131 1209 62 1600 43 Data for Problem 1 Graph the data points Find a third order Lagrange polynomial u f3p suitable for interpolation over the entire range of densities Plot the polynomial on the same graph as the data points Calculate the true error fp f3p where u p is the true function at the two remaining data points Estimate the average speed of traffic when the mile and a half long tunnel has 75 cars uniformly spaced in it Suppose the speed limit in the tunnel is 35 mph Estimate the uniform traffic density throughout the tunnel at which traffic tends to ow at the speed limit The traffic ow q vehicleshr through the tunnel under steadystate conditions uniform density and constant speed is the product of traffic density p vehiclesmile and vehicle speed u mileshr Generate a graph of q vs p for the tunnel using the Lagrange polynomial from Part b Traffic engineers refer to this as the Fundamental Diagram of Traffic Flow g The maximum traf c ow is called the capacity of the road Estimate the tunnel s capacity What should the posted speed limit be to achieve this capacity 2 The eld strength of the earth39s gravitational eld g diminishes with distance h from the surface of the earth The following table shows the eld strength at several distances Distance h Field Strength g Distance h Field Strength g 10 m Nkg 10 m Nkg 0 981 20 058 25 507 30 030 5 309 40 019 10 149 50 012 Data for Problem 2 a Plot the data points b An interpolating polynomial is required to estimate g from zero to 20 gtlt106 In Find the Lagrange interpolating polynomial f7h which passes through all the data points and the Lagrange interpolating polynomial at which passes through the first ve data points Plot both over the interval 0 S h S 20gtlt106 m and comment on which is more suitable for interpolation c Estimate the distance from earth where a 1 kg mass weighs 2 Newtons ie where the gravitational eld strength is 2 Nkg using the two interpolating polynomials form Part b d The true function relating g and h is Gm g f h EZ RE h G 67X103911Nrm2kg2 mE 60gtlt1024 kg RE 64gtlt106 m where gravitational constant mass of the earth radius of the earth Plot the true function h on the same graph as the data points and the interpolating polynomials at and f7h over the interval 0 S h S 50gtlt106 m Comment on the suitability of either interpolating polynomial over the entire interval 0 S h S 50gtlt106 m e Calculate the true errors in the two estimates of g in Part c ry to nd an tnterpolaung polynomlal that aeeurate estrrnates of g over data potnts from the true funeuon g Use an approprtate Lagrange r from the earth s surf aee where y the enttre tnteryal 0 lt g M tf nterpolattng polynomlal to est ou Welgh nehal ean be used to provlde reasonably h 50x105 rn Generat more you thtnk thts wlll help trnate the dstanee fof your wetght on earth h Caleulate the errortn the esurnated dtstanee tn Partg wtth a bat speed of 70 rnp eorrespond to the arnount th w en the underswtng ls zltnehes ls ll Souree Falk SportFaets F w Roth The euryes shown below are the txajectones of a batted baseball thrown at 85 rnph h dreeted upward by 10 de r the bat is swung undert The vanous curves g es he ball For example the ease er lustratedtn the upper eorn 1991 w 20 in 15 in E E z E 10 tn 5 25 n 035 ln 1 so 05 in a tn 350 n 1 2m DlsuNcE N FEEY Graphs for Problem 3 a Seleet one of the euryes and gener from th ate a set of ughtly spaeed data potnts sarnpled e eurye Represent the aetual eurye by a pteeewt e lmear approxlmanon through the data potnts Grap e apprornrnaung funetton whteh should appear srnooth as a dotted eurye wtthout showtng the data potnts Label rt HjD b For the eleet Part a generat sarne s ed tn potnts no rnore than 10 that eaptu e a re uee res the baste proftle ofthe eu d set of sarnpled data rye e Try to flt a Lagrange tnterpolattng polynomlal D over the enttre tnteryal 4 44 mlmml w e e polynomlal D andthe data potnts on the sarne gaph 4 An electronics lab experiment was designed to compare the measured voltagecurrent relationship for a light emitting diode with calculated values based on the known nonlinear resistance characteristic of the diode Results are given in the table below Source Monaghan Computers in Education JanMar 1998 VOLTAGE VS CURRENT FOR A LIGHT EMITTING DIODE LAB MEASUREMENTS CALCULATED VALUES VOLTAGE CURRENT VOLTAGE CURRENT volts milliamperes volts milliamperes 1421 0095 1420 0118 1441 0155 1440 0186 1446 0281 1460 0295 1487 0461 1480 0468 1512 0808 1500 0741 1535 1390 1520 1173 1556 2350 1540 1859 1557 4080 1560 2944 1596 6720 1580 4664 U V Table for Problem 4 Plot the measured and calculated current vs voltage data points on separate graphs Find second third and fourth order Lagrange interpolating polynomials for each set based on the following data points MEASURED ORDER DATA POINTS 2 142100951512 080815966720 3 142100951487 04611535139015966720 4 142100951446 02811512 0808 1556 235015966720 CALCULATED ORDER DATA POINTS 2 14200118 1500 0741 15804664 3 1420011814800 4681520117315804 664 4 1420011814600 2951 500 0 7411540185915804 664 e P1ot the po1ynornra1s on the sarne graphs wrth the data pornts Comment on the resu1ts d Denote x v as the funcuon re1aung the aetua1 rneasured eurrent and vo1tage Ca1eu1ate the average error for eaeh rnterpo1aung po1ynornra1v 1 e r Em 1 72V mm n 234 3 Repeat Part d wrth x v as the funeuon re1atrng the aetua1 ea1eu1ated eurrent and vo1tage gure 1800220 H 150 200 250 300 350 9 Cam Ang e Degvees u an 100 Frgure for Problem 5 a Generate two sets of data pornts 9 x1 one for earn ang1es a between 0 and 180 d ree the h f e 1es between 220 and 360 degree Ch e pornts so at L range rnterpo1atrng po1ynornra1s ean be used to approxxmate both pro les P1otboth sets of data pornts on the sarne graph 05 b c d Find the lowest order Lagrange interpolating polynomials that describe both profiles reasonably well They do not have to be the same order Plot the polynomials from Part b on the same graph as the data points Select several additional data points from both pro les and calculate the true error f9 fn9 where n is the order of the polynomial at each of those points A car buyer is nancing the purchase of a 35000 automobile The bank has several loan packages to choose from The choices are summarized below Loan Period Annualized Interest rate Monthly Payment 71 months 139 A 24 50 153550 36 55 105686 48 60 82186 60 65 68482 72 70 59672 Table for Problem 6 Plot the data points nkAk k 12345 Find the 43911 order Lagrange interpolating polynomial A f4n through the data points and plot it on the same graph The bank is willing to write the loan for any number of months between 24 and 72 Estimate the monthly payment if the car buyer chooses a three and a half year loan The interest rate on the loan increases with the duration of the loan Use the table to interpolate 13942 the annualized per cent interest rate for a fortytwo month loan The formula for calculating the monthly payment A given the amount nanced P the annual per cent interest rate 139 and the loan duration 71 in months is AfPz39 nP 15W115Wquot 1732100 1 Find the true error in the estimated monthly payment 42 calculated in Part c ie 35000 i42 42 42 The car buyer can afford monthly payments of 770 Use the graph A f4n to estimate to the nearest month the loan duration requiring this amount of payment g Find the true loan duration to the nearest month when the monthly payment is 770 and compare the result to the answer in Part f Hint You will need to solve A P 1 n for n by trial and error since 1 is unknown as well 7 The estimated net lifetime bene ts in in ation adjusted 1993 dollars for social security recipients born in different years is tabulated below Source US News and World Report April 20 1998 Birth Year Net Lifetime Benefits 1885 29000 1900 51000 1915 60000 1930 14000 1945 18000 1960 33000 1975 32000 1985 31000 Data for Problem 7 a Graph the data points b Find the seventh order Lagrange interpolating polynomial passing through the data points and plot it on the same graph c Estimate the year of birth for a social security recipient to have net benefits of zero d Estimate the maximum net lifetime social security benefits and the year of birth for recipients receiving the maximum 8 Consider the 3 points from the function fx ex in the table below 1 x y fx 0 1n1 1 1 In 12 12 2 1n7 7 Table for Problem 8 a Find the first order Lagrange polynomial f1x to be used for estimating fln5 0 0 b Find the second order Lagrange polynomial f2x to be used for estimating fln5 c Plot the given data points and both interpolating polynomials on the same graph d Calculate the interpolation errors Riln5 fln5 x 139 1 2 Railroad freight traffic measured in billions of tonmiles since 1965 is presented in tabular form Source Assoc of American Railroads Year 1965 1970 1975 1980 1985 1990 1995 lTrafficl 700 750 750 925 850 980 1275 Data Points for Problem 9 Estimate the year between 1975 and 1985 when the freight traffic was a maximum and the amount of traffic in that year The following data points were obtained from a function y x x 0 1 2 3 4 5 6 7 8 9 10 fx 10 8 6 4 2 0 2 4 6 8 10 a Plot the points b Find the lowest order Lagrange interpolating polynomial through the entire set of data points and plot it on the same graph c Suppose the middle data point was incorrectly obtained as 5 y where y 05 Find the lowest order interpolating polynomial x n S 10 through the entire set of data points and plot it on the same graph The interpolating polynomial may be expressed in standard form d Find the second order Lagrange interpolating polynomial fzx through the end points and the new middle data point and plot it on the same graph e Assuming the function fx is linear calculate the true error at x 04 and x 475 resulting from the use of both polynomials f Repeat Parts c d and e fory 1 and 2 g Comment on the results ARBITRARY UNKNOWNS The echelon form of the augmented matrix con rms the existence of arbitrary unknowns ie a consistent system of equations in which one or more variables can be chosen arbitrarily There are several ways to establish if indeed a certain variable can be included in the subset of arbitrary unknowns A few simple examples illustrate the point Example 1 For the system of equations below establish the existence of l arbitrary unknown and determine if x y and 2 can each be arbitrary x y z 0 2x y z 3 x 4y 42 3 x 2y 22 3 1 1 10 1 1 10 1110 I I 2 1 13 0 3 33 011 1 Ab I N I N I 1 4 4 3 0 3 3 3 0 0 0 0 1 2 2i3 0 3 3i3 oooio The 4 by 4 system has been reduced to a 2 by 3 system That is there are only 2 linearly independent equations in the 3 variables unknowns and hence one of the variables is arbitrary These equations are xyz0 Now suppose we try to make 2 the arbitrary unknown This requires that solutions for x and y be expressed in terms of z This can be attempted in one of two ways One approach is to simply transfer all terms involving 2 over to the right side of the equations and consider 2 as a parameter This yields xy z y l z All that remains is substituting y l 2 into the rst equation and then solving for x The nal solution can be expressed as x 1 y lz z arbitrary Clearly an in nite number of solutions exist since there is a different solution for each arbitrarily assigned value for z A slightly different approach involves reformulating the reduced equations from the echelon form as y A A A A l l A z Ax b where A b 0 l l z ie as a new system in matrix form with modi ed coef cient matrix A constant vector 13 and vector of unknowns i The identical solution as given above is easily obtained from the modi ed augmented matrix M 3 154 3 354 xl ylz zarbitrary IN The second approach is somewhat more instructive because the modi ed coef cient matrix A determines whether the variables placed on the right hand side 2 in this example are arbitrary When A which will always be a square matrix is nonsingular has a unique solution and the variables moved to the right hand side are indeed arbitrary Consider what happens when x is selected to be the arbitrary variable The reduced 2 by 3 system with x as the arbitrary unknown becomes and attempting to solve for a solution by GaussJordan gives M 1 iii N 3 SEQ In this case a unique solution for y and z in terms ofx was not possible This comes as no surprise because A is clearly singular This con rms what was already known from the previous case where 2 was arbitrary but x was constrained to be 1 In fact the GaussJordan solution above implies that x l 0 for the equations to be consistent For larger systems with arbitrary unknowns it is essential to check the modi ed coef cient matrix A before continuing on with a GaussJordan or backward substitution under the assumption that one or more speci c variables can be arbitrary The following system demonstrates this point x1 x2 x3 x 4 x5 2x6 3 x2 x4 3x6 0 x1 x3 x4 2x5 4x6 4 2x1 x2 x3 4x6 1 x1 3x2 2x3 x4 2x6 1 x2 x3 2x4 x5 8x6 0 11111 23 01010 30 001 3 11150 Aline 0001 1 3 1 00001051 0000000 The original 6 by 6 system of equations has been converted to a new 5 by 6 system with a single arbitrary unknown From the echelon form these equations are x1 x2 x 4 3x6 0 x3 3x4 x5 1 1x6 0 x4 x5 3x6 l x 1 Clearly x5 is not arbitrary The modi ed matrix A obtained from the columns of the echelon form with the x5 column omitted is shown below Quite obviously its singular as expected confirming that x5 is not arbitrary 1111 2 0101 3 13001 311 0001 3 00000 However there may be other nonarbitrary variables in addition to x5 which are not as obvious For example suppose we proceed to solve the reduced 5 by 6 system using GaussJordan with x2 selected as the arbitrary unknown Observe what happens 1 1 1 1 2I3 x2 1 1 1 1 2 3 x2 00 1 0 3i x2 01 3 1 11 0 ME 0 1 3 1 11E 0 0 0 1 0 3 x2 00 1 1 3i 1 00 1 1 3 1 0 0 0 1 0l 1 0 0 0 1 0 1 1 1 1 1 2 3 x2 1 1 1 1 2 3 x2 01 3 111i 0 01 3 111 0 0 0 1 0 3i x2 0 0 1 0 3 x2 0 0 0 1 0i 1x2 0 0 0 1 01 x2 0 0 0 1 0 1 0 0 0 0 0 x2 The bottom row of zeros in the rst 5 columns signi es that a solution for x5 is impossible when x is chosen to be arbitrary and the GaussJordan method terminates without a solution Furthermore for consistency the last row implies that x must be zero further eVidence it can t be arbitrary The 41 row represents the equation x5 l x2 which is consistent with the last row of the echelon form which states that x5 1 In problems of this type the prudent thing to do is verify that the modi ed coef cient matrix A is nonsingular before proceeding to nd a solution In the preVious example when x was assumed to be arbitrary A became 1111 2 0010 3 A01 3 111 001 1 3 00010 From MATLAB its easy to verify that A is singular and therefore x2 should not be chosen as arbitrary A l l l l 2 0 0 l 0 3 0 l 3 l l l 0 0 l l 3 0 0 0 l 0 EDU detA ans 0 The same approach applies when more than one variable is arbitrary To illustrate consider the system of equations x1 x2 x3 x4 3x5 3 x1 x2 x3 3x4 x5 1 2x2 x3 x 4 2x5 4 x1 2x3 5x5 1 2x1 x2 x3 3x 4 4x5 0 3x1 2x2 6x4 3x5 1 The echelon form of the augmented matrix has three rows consisting entirely of zeros The rst three nonzero rows are 111 1 33 l 0111 22 0011 2io Ab 1 1 1 1 3 3 1 1 1 3 1 1 0 2 1 1 2 1 0 2 0 5 1 2 1 1 3 4 0 3 2 0 6 3 1 EDU rrefAb ans 1 0 0 2 1 1 0 1 0 0 0 2 0 0 1 1 2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 In the 3 by 5 system of equations that correspond to either echelon form there must be 2 arbitrary unknowns To check if say x4 and x5 can be arbitrary we look at the modi ed coef cient matrix A that results when the columns for x4 and x5 are removed from the rst echelon form b H COD I Ol D I 1 1 1 V Since A is a nonsingular matrix there is a unique solution to Ag h 3x4 3x5 where 3 x2 and b 2 x42x5 x4 2x5 1 l 0 3x43x5 l2x4x5 The solution is g All 0 l l 2 x4 2x5 2 0 0 1 x4 2x5 x4 2x5 ie x1 1 2x4 x5 x2 2 x3 x4 2x5 x4 arbitrary x5 arbitrary Since we know from the previous solution that x is not arbitrary it should come as no surprise that any 3 by 3 submatrix formed from the rst ve columns of either echelon matrix minus the zero rows is destined to be singular if it excludes the second column the one corresponding to x2 The resulting 3 by 3 matrices obtained from the rst echelon form are given below and the reader should verify that they are all singular 1 1 3 1 3 1 1 2 x1 and x2 columns removed 0 1 2 x2 and x3 columns removed 1 1 2 0 1 2 l l 3 l l l 0 1 2 x2 and x4 colums removed 0 1 1 x2 and x5 columns removed 0 l 2 0 l The row reduced echelon form is even more explicit as to why x2 cannot be arbitrary It is clear from this echelon form that x 2 and hence not arbitrary Furthermore removing the x2 column from the row reduced echelon form leaves the following matrix 10 2 11 I 00 0 01 2l0 Regardless of which additional column is removed for the 2quotd arbitrary variable the resulting 3 by 3 modi ed coef cient matrix A will be singular con rming that x cannot be one of the two arbitrary unknowns In summary either echelon form of the original augmented matrix will reveal the existence of arbitrary unknowns The original m by n system of m equations in n unknowns will be reduced to an ml by n system where m1 S m indicating that m m1 equations from the original system were redundant Ifml is less than n the existence of n m1 arbitrary unknowns is assured A particular subset of n m1 unknowns is arbitrary provided the ml by m1 submatrix of the echelon matrix obtained by removing the columns corresponding to the n m1 unknowns is nonsingular The row reduced echelon form in MATLAB will identify the arbitrary variables directly ie any row with a single 1 in the rst n columns implies the variable associated with that column is not arbitrary For example in the following 4 by 6 rowreduced echelon form x1 and x5 cannot be arbitrary unknowns 10000 010534E1 rrefAl 00127155 000010E8 000000 Chapter 18 Polynomials 181 Roots MATLAB performs numerous operations involving polynomials The most common situation dealing with polynomials is nding its roots A polynomial is uniquely de ned by its coefficients MATLAB represents a polynomial by a vector containing the coef cients of the powers in descending order For example the polynomial 71 1706 61xquot aHxquot ch a0 is saved in MATLAB as a vector p where p anaH a1a0 If a coef cient is zero it must appear in the vector A er the coef cient vector is de ned the roots of the polynomial are returned in a column vector from a call to roo t s p The roots of a polynomial are found in the example below The polynomial is evaluated at the roots where it should be zero and then plotted to verify it passes through the real roots Example 1811 p110429865 Define coefficient vector for polynomial px x4 10x3 42x2 98x 65 prootsrootsp Find roots of px xproots Set x equal to row vector of roots of px0 yplxA4 p2xA3 p3xA2 p4x p5 Evaluate px at each of its 4 roots xilinspace010250 Create xi vector for evaluating yipxi yiplxiA4 p2xiA3 p3xiA2 p4xi p5 yipxi plotxiyi Plot ypx vs x hold on axis0 6 50 50 Set scale for plot plot0 100 0 k Plot xaxis xlabel x ylabel y title ypx vs x p l 10 42 98 65 20000 30000i 20000 0000i 50000 10000 A 109012 00568 00568i 00568 00568i 03411 0 ypx vs X The reveise process ie obtaining the polynomial coef cients from the roots is easily accomplished using the function 39p o l y Example 1812 pcoeffspolyproots Find coeff s of polynomial with roots proots pcoeffs 10000 100000 420000 980000 650000 182 Multiplication Multiplying two polynomials in MATLAB is as simple as de ning the coef cient vect01s and then performing a convolution using the 39c o nv function Example 1821 p12 54 First polynomial is px x4 2x3 5x 4 ql 0 2 0 Second polynomial is qx x3 2x uconv P I q p 1 2 5 4 q 1 o 2 o u 1 2 7 8 10 8 o Executing Example 1821 shows that pxqx x4 2x3 Sx 4x3 2x x7 2x6 7x5 8x4 10x38x The roots of the product polynomial and the individual polynomials must agree This is easily ven39 ed Example 1822 rprootsp Find roots of px rqrootsq Find roots of qx rurootsu Find roots of uxpxqx 31774 18558 O6784 14142 14142 31774 14142 18558 14142 O6784 183 Addition Addition of polynomials is a bit more complicated than multiplication when the polynomials involved are not the same order That is MATLAB cannot add polynomials of diiTerent order directly because the coef cient vectors representing them are not the same size The solution is to pad the lower order polynomial coef cient vector with leading zeros until the two coef cient vectors are the same size The function 39polyadd in Example 1831 below adds two polynomials Several test cases are included after the function Example 1831 function hpolyaddfg This function adds two polynomials If the polynomials are not the same order it pads the coefficient vector of the lower order polynomial with zeros so the two vectors can be added if lengthflengthg Check if polynomials are the same order hfg elseif lengthfgtlengthg Check if f is higher order than g difflengthflengthg ggzeros1diff g Pad g with diff zeros hfgg Compute sum else g is higher order than f difflengthglenghtf ffzeros1diff f Pad f with diff zeros hffg Compute sum end if lengthflengthg disp39 disp Case 10 gl 2 3 4 Define gx x3 2x2 3x 4 h1 2 3 4 Define hx x3 2x2 3x 4 zpolyaddgh Add polynomials g and h disp39 disp Case 239 f1 3 2 o 5 Define fx x4 3x3 2x2 5x g4 1 6 Define gx 4x2 ypolyaddfg Add polynomials f and g disp disp Case 3W a1 ax1 b10001 bxx41 ypolyaddab Case 1 g l 2 3 4 h l 2 3 4 z O O O 0 Case 2 f l 3 2 O 5 g 4 l 6 y l 3 2 l 11 Case 3 a l b l O O O l y 1 o o o 0 Note the use of the local variable blif f 39 in the lnction workspace Despite the fact a builtin MATLAB lnction named iii f f 39 exists there is no problem assigning it a value and refering to it in the next statement Recall that MATLAB rst looks in the current worlspace for named valiables before checking if its a builtin Jnction 184 Division Division of polynomials is achieved by the deconvolution function deconv In the example below the polynomial N x 4x5 2x3 x2 x l is divided by the lower order polynomial Dx x2 2 The output is returned in two vectors The fust is the coef cient vector of the quotient and the second produces the remainder polynomial coef cients In general Noe Roe Doc Q Doc where Qx is retumed as zero if the order of N x is less than the order of Dx Example 1841 N4 o 2 1 1 1 Define Nx 4x5 2x3 x2 x 1 D1 o 2 Define Dx x2 2 AqArdeconvND Find NxDx N 4 o 2 1 1 1 D 1 o 2 Aq 4 o 10 1 Ar o o o o 21 1 A er running Example 1841 the following result is apparent 4x52x x2xl 4x310x1 zlx l x 2 x 2 The above result is easily veri ed in MATLAB by showing that NOC Q06 Dx Rx Example 1842 N Display Nx coefficients NcheckpolyaddconvAqD Ar Check result from Example 1841 N 4 O 2 l l l Nicheck 4 O 2 l l l The next example demonstrates the case when the numerator polynomial is lower order than the denominater polynomial Example 1843 A1 2 3 Define Ax x2 2x 3 B1 o 1 o 51 Define Bx x4 x2 5 CqCrdeconvAB Find AxBx A l 2 3 B l O l O 5 Cq 0 Cr 1 2 3 185 Differentiation Differentiation of polynomial functions is a straightforward procedure one that MATLAB implements with the function po 1 yde r Example 1851 f13 Coefficient of highest order term Coefficient of next highest order term 1 f52 Coefficient of zero order term ie the constant ffx 3x45x3x2 dfdxpolyderf Find dfdx coefficient vector f 3 5 O l 2 C1de 12 15 0 1 Example 1852 f1 0 3 2 fx x3 3x 2 gl 4 5 gx x2 4x 5 fgconvfg Find fxgx coefficient vector dfgdxpolyderfg Differentiate the product fxgx u1convfpolyderg Find fdgdx u2convgpolyderf Find gdfdx dfgdxpolyaddu1u2 Differentiate fxgx using product rule f 1 o 3 2 g 1 4 5 fg 1 4 8 1o 7 10 difgidx 5 16 24 2o 7 difgidx 5 16 24 2o 7 186 Evaluation Polynomials are usually evaluated at one or more discrete points over some interval of interest In MATLAB 39polyval is the builtin function for doing this To illustrate polynimial evaluation consider an annual deposit A which earns interest at 139 per annum The lture worth of the account a er 71 years is F AA1z39A1z 2 A1z H A1z39quot 1 A1 iquot Auquot uquot 1 uquot392 u2 u 1 where u li The following example computes the account balance a er 71 years by evaluating the polynomial function above Example 1861 A1200 Annual deposit i10 Annual interest rate per cent ii100 Convert i to decimal value u1i n15 Number of years pones1n1 Coefficient vector of polynomial FApolyvalpu Compute future amount A 1200 i F 43140e004 The 39polyval39 function provides a quick way to establish the values of a polynomial so it can be plotted Example 1862 A1200 ilinspace025500 Interval of i values for plot ii100 Convert array of i values to decimal numbers u1i n15 Number of years pones1n1 Coefficient vector of polynomial FApolyvalpu10000 Future worth10000 i100i Convert i back to per cent plotiF Plot E vs i title Future Worth of 16 Annual 1200 Deposits Vs Annual Interest Rate i0 xlabel Interest Rate per Annum ylabel Future Worth x 10000 Future Worth of 16 Annual 1200 Deposits Vs Annual Interest Rate i Future Worth X 10000 I I 0 5 10 15 20 25 Interest Rate per Annum Example 1863 aconv1 0 11 2 Coefficient vector of fxx2lx2 frootsrootsa xi20013 Values of x to evaluate fx at yipolyvalaxi yifxi plotxiyi Plot yfx vs x hold on P10t2300 k Plot xaxis xlabel x ylabel y title fx vs x a l 2 l 2 firoots 20000 10000 lOOOO fx vs x 187 Rational Polynomials There are times when the ratio of two rational polynomials needs to be expressed x3 1n part1al fraction form For example the express1on can be x 5x 7x2 expanded into the sum of 3 terms using the MATLAB re s idu e function Example 1871 n1 3 Numerator polynomial nx x 3 dl 5 7 2 Denominator polynomial dx x3 5x2 7x 2 rpkresiduend Find residues and poles of nxdx n l 3 d l 5 7 2 r 02764 lOOOO 07236 p 26180 20000 O3820 k H From the results obtained in the above example it follows that x3 02764 1 I 07236 x35x27x2 x26180 x2 39 x0382 When the order of the numerator polynomial is greater than or equal to the order of the denominater polynomial the last output k will not be an empty vector Instead it will contain the coef cients de ning the non fractional term in the expansion Example 1872 nld dln Reverse nx and dx rlp1klresiduenld1 Find partial fraction expansion n1 l 5 7 2 d1 l 3 r1 1 pl 3 kl l 2 l x35x27x2 1x12x1 x3 x3 In this case The 39residue function can also be used to sum the individual terms of the partial fraction expansion There must be 3 input vector arguments r p and k an empty vector in the absence of a nonfractional term Example 1873 IIP ndresiduer p Combine partial fraction terms into nd r 02764 l0000 07236 26180 20000 03820 n 00000 10000 30000 d 10000 50000 70000 20000 r1plkl n1dl residuerl p1 k1 Combine rlpl k1 polynomial r1 1 pl 3 kl l 2 1 n1 l 5 7 2 d1 l 3 188 Curve Fitting There are a variety of ways to represent a set of empirical data One method Least Squares Regression is based on the use of low order polynomials to t the data in the least squares sense A least squares t is a unique curve with the property that the sum of the squares of the differences between the tted curve and the given data at the data points is a minimum MATLAB generates least squares polynomials for a set of data by using the function bo1yfit It accepts vectors of x and y data and a scalar n as the order of the polynomial approximation The target heart rate during exercise for diiTerent ages is tabulated below IXge 20 25 30 35 40 45 50 55 60 65 70 nga HemRate 150 146 142 139 135 131 127 124 120 116 112 The data is plotted and a rst order least square polynomial is found and presented on the same graph in the next example Example 1881 x20z5z70 Create vector of age values y150 146 142 139 135 131 127 124 120 116 112 Create vector of target heart rate values plotxy Plot raw data a1polyfitxy1 Find coefficients of least square line xilinspace2070250 Create vector of x values for plotting yipolyvala1xi Evaluate points on least square line for plotting hold on v20 70 100 150 plotxiyi axisv xlabel Age years ylabel Target Heart Rate x 20 25 3O 35 4O 45 5O 55 6O 65 70 y 150 146 142 139 135 131 127 124 120 116 112 a1 O7527 1649636 Target Heart Rate I 20 25 30 35 40 45 50 55 60 65 70 Age years A higher order polynomial curve is warranted in the next example where air pressure and the speed of sound valiations with altitude in the troposphere are tabulated Altitude h it Pressure p psi Speed of Sound v sec 0 147 11164 5000 122 10971 10000 101 10774 20000 676 10369 30000 437 99485 36000 331 96875 Example 1882 h0 5 10 20 30 36 Altitude data pl47 122 101 676 437 331 Pressure data v11164 10971 10774 10369 99485 96875 Speed data p2polyfithp2 Find quadratic fit coefficients for p hilinspace040500 Create h data for evaluation of 2nd order fit pipolyvalp2hi Compute quadratic fit p values v2polyfithv2 Find quadratic fit coefficients for v vipolyvalv2hi Compute quadratic fit v values subplotl2l Divide figure Window for 2 side by side graphs and make the left side active hold on plothp Plot p vs h table data in left side plothipi Plot quadratic fit hi vs pi in left side xlabel Altitude h ft x 1000 ylabel Air Pressure p psi title Air Pressure Vs Altitude subplot122 Make right side active hold on plothv Plot v vs h table data in right side plothivi Plot quadratic fit vi vs hi in right side xlabel Altitude h ft x 1000 ylabel Speed of Sound v ftsec title Speed of Sound Vs Altitude Air Pressure Vs Altitude Speed of Sound Vs Altitude 16 1120 a 1100 1080 1060 v ftsec 1040 1020 Air Pressure p psi 1000 Speed of Sound 2 l l 940 l l l l 0 10 20 30 40 0 10 20 30 40 Altitude h ft x 1000 Altitude h ft x 1000 Chapter 14 Function M les 141 M File Construction Rules There is another type of mFile besides a MATLAB script le It is a MATLAB function mFile Typically a MATLAB function le contains a sequence of MATLAB commands which process speci c information passed to it in the form of input variables and generate output based on the information A function mFile behaves like a black box from the user s standpoint because all that is evident are the variables passed to it from the MATLAB Workspace and those which are returned from it to the Workspace Intermediate variables created in the function mFile are not stored in the MATLAB Workspace and hence will not overwrite any variables of the same name which happen to be present in the Workspace Function mFiles are created with the MATLAB EditorDebugger and saved as text les with a m extension similar to script les All MATLAB built in mathematical functions are function m Files A simple function mFile called f loan m is shown in Example 1411 Its main purpose is to calculate amathematical function value based on values passed to it from a script le called Loanilnterest m and return the computed value to the calling program The calling script le Loanilnterest m is included as well in order to see where the function calls occur Example 1411 function iR floanAAPPNNii floan is a function called by the script file LoanInterestm It accepts a monthly payment AA loan amount PP a loan duration NN in years and an interest rate ii It returns a single value determined from evaluation of the function The script file Loanlnterestm is a root solving program which attempts to find the root iR of fAPNi0 The user selects from one of two root solvers the Bisection Method or the Method of False Position The root iR found in LoanInterestm is the annualized interest in percent dePdePdePdePdP z monthsl2NN Convert from years to months if ii0 Check if ii not equal 0 xlii iRAAPPiixANmonthsxANmonthsl Calculate function value when ii not equal 0 else ii equals 0 iRAAPP 1Nmonths Calculate function value when ii equals 0 end if ii0 Check if ii not equal 0 LoanInterestm A script file that calculates the iterest rate on a loan It calls function floanm with inputs PAN and i PzLoan amount AzMonthly installment NzDuration of loan expressed in years izAn interest rate decimal value iR Annual per cent interest rate on loan dePdePdePdeP o clearclc P35000 Loan amount A750 Monthly payment 5 Loan duration in years method F iLinitial0 iUinitial1 fiLfloanAPNiL fiUfloanAPNiU count0 i eAeSTOP1 While countlt100ampeAgteSTOP ioldi Calculate next i and fi if method B iiLiU2 iiUfiUiLiUfiLfiU end if fifloanAPNi inewi Determine new endpoint and fi at new endpoint if fiLfilt0 iUi fiUfloanAPNiU else iLi fiLfloanAPNiL end if Caclulate approximate relative error and increment counter if countgt0 eAabsineWioldinew end if countcount1 end While iR1200i Calculate annual per cent interest rate fprintf nn fprintf Loan Amount 80f nMonthly Payment 71f nLoan Period 40f years n PAN fprintf Annual Interest Rate 62f per cent n iR Enter Loan Amount 40000 Enter Monthly Installment 900 Enter Loan Duration in Years 5 Enter quotBquot for Bisection Method or quotFquot for False Position F Loan Amount 40000 Monthly Payment 9000 Loan Period 5 years Annual Interest Rate 7 1250 per cent Another Problem Y or N N Several important characteristics of function m Files are The rst line must contain the word function followed by its calling syntax see function f loan m The input and output variables on the function declaration line are local to the function meaning they are not accessible to the calling le which may be a script le or possibly another function le They are merely temporary variables with values inside the function39s workspace In the previous example input variables AA PP NN and ii assume values passed to the function from the calling script le and the value of output variable iiR computed inside the function is passed back to the calling script le Input variables can be modi ed inside the function ie appear on the left side of an assignment statement but the new values cannot be retumed to the calling le N Function names have similar constraints to MATLAB vanable names E The le named used to save the function mFile should be identical to the name appearing in the function declaration line although its not required MATLAB references the le name stored on disk when executing a Jnction 4 A function m File terminates after executing the last statement or whenever a r e t u r n command is encountered 5 Function mFiles can call script les and other functions which are referred to as subfunctions The subfunctions are comtructed in the same was as a function mFile and are appended to the end of the primary function in any order The next example shows a sequence of commands entered in the Command window A function m le fl m is called with arguments A B and i which have been assigned numerical values The function 39f l m shown below contains a reference to a subfunction mFile nois e m The 39noi se m subfunction declaration line appears at the end of the lnction 39f l m listing Example 1412 uuflABi end i 2 09681 y 39681 u 39681 2 06199 y 56199 u 95880 2 06408 y 76408 u 172287 function yf1ABx Linear function with noise component yABxnoiseA function znoiseC Noise function zCrand The subfunction hoise m is invoked by the function fl m which passes it the value of variable A namely 2 to be substituted for variable C wherever it appears in noise m Hence the rst output variable calculated is the subfunction variable 2 whose value is passed back to the calling function fl m as the noise component in the calculation of y Finally the numerical value of y is returned to the Command window program to update variable u The comment lines following the function declaration deserve careful consideration They will be displayed whenever a help 39 function name command is issued usually from the Command window For example to obtain a description of functions 39floan39 and fl 39 see Example 1413 Example 1413 help floan Get description of function floan floan is a function called by the script file 39Loanilnterestm39 lt accepts a loan amount P a monthly payment A a loan duration N in years and an interest rate i It returns a single value determined from evaluation of the function The script file 39Loanilnterestm39 is a root solving program which attempts to find the root iR of fAPNi0 The user selects from one of two root solvers the Bisection Method or the Method of False Position The root iR found in 39Loanilnterestm39 is the annualized interest in percent 142 Input and Output Arguments There are a number of rules governing the use of input and output arguments in MATLAB function calls Several important features are listed below Nt I E 4 A function mFile can have zero input and zero output arguments A function can be called with fewer input and output arguments than appear in the function declaration line The number of input and output arguments passed to the mFile function can be determined inside the function by using the MATLAB functions nargin and n a rgout39 respectively A function with outputs in the function declaration line is not required to return values for all of its outputs Any output not assigned a value in the body of the m File function will not retum a value back to where its called from Example 1421 function x0yOoriginclosexly1x2y2x3y3x4y4 This functions accepts up to 4 input points xiyi and finds the one closest to the origin The coordinates of the closest point to the origin are returned in the outputs x0 and y0 switch nargin case 0 disp No points searched case 2 x0xly0yl case 4 d1sqrtxlxl ylyl d2sqrtx2x2 y2y2 if d1 lt d2 x0xly0yl else x0x2y0y2 end if dlltd2 case 6 d1sqrtxlxl ylyl d2sqrtx2x2 y2y2 d3sqrtx3x3 y3y3 if d1mind1d2d3 x0xly0yl elseif d2mind1d2d3 x0x2y0y2 else x0x3y0y3 end if d1mind1d2d3 case 8 dlsqrtxlxl d2sqrtx2x2 y2y2 d3sqrtx3x3 y3y3 d4sqrtx4x4 y4y4 if d1mind1d2d3d4 y1y1 x0xly0y1 mind1d2d3d4 Y r elseif d3mind1d2d3d4 x0x3y0y3 else x0x4y0y4 end if d1mind1d2d3d4 end switch nargin The following MATLAB commands create a set of 4 data points and then calls the function 0 r i g i mic l o s e m 39 to nd the closest point to the origin x14randl4y14rand14 xlx1 y1y1 x2x2 y2y2 x3x3 y3y3 x4x4 y4y4 xnearynearoriginclosex1ylx2y2x3y3x4y4 x 00136 07895 08996 00022 y 05163 06244 02652 03990 xnear 00022 ynear 03990 The function in Example 1422 has no outputs It results in the drawing of a green rectangle of specified size and location determined by the numerical values passed to its inputs The inputs fix the lower and upper comer vertices subject to a set of constraints A script le that calls the function several times is included Example 1422 function roadmedianxAyAxCyC MATLAB function to draw medians xAyA and xCyC are the coordinates of the lower left and upper right corners OixAxC31000 and OiyAyC31000 if xAgt0ampyAgt0ampxClt1000ampyClt1000 Check if vertices are within limits hold on v0 1000 0 1000 Set scale xBxAyByC Set coordinates of pt B upper left xDxCyDyA Set coordinates of pt D lower right xxA xB xC xD yyA yB yc yD fillxy g Draw green rectangle axisv else warning 0ne or More Vertex Coordinates Outside of LimitsU end if xAgt0ampyAgt0ampxClt1000ampyClt1000 xll0yll0 x1u100y1u150 x2l350y2l200 x2u500 y2u500 x3l750 y3l400 x3u900y3u800 roadmedian xll yll xlu ylu roadmedian x21 y2l x2u y2u roadmedian x31 y3l x3u y3u roadmedian x31 y3l x3u y3u300 Warning One or More Vertex Coordinates Outside of Limits gt In CMATLABRllW0rkMatlabicourseroadimedianm at line 20 1000 I I I I I I I I I 0 100 200 300 400 500 600 700 800 900 1000 The next example contains a function which draws a traf c signal face at a location set in the calling script le The script le also controls the face colors which are lit by Virtue of input argument values passed to the function The function is listed below followed by the script le Example 1423 function TrafficSignal xyvisiblegvisibleyvisibler Traffic Signal plots traffic signal with green yellow or red visible at location xy BEGIN PLACEMENT OF TRAFFIC SIGNALS xtsgx 10 x coordinate of traffic signal green light ytsgy 10 y coordinate of traffic signal green light xtsyxtsg 2 x coordinate of traffic signal yellow light ytsyytsg y coordinate of traffic signal yellow light xtsrxtsy 2 x coordinate of traffic signal red light ytsrytsy y coordinate of traffic signal red light xtsAxtsg 2 x coordinate of traffic signal lower left corner ytsAytsg 8 y coordinate of traffic signal lower left corner xtsBxtsr 2 x coordinate of traffic signal lower right corner ytsBytsA y coordinate of traffic signal lower right corner xtsCxtsB x coordinate of traffic signal upper right corner ytspytsB 16 y coordinate of traffic signal upper right corner xtsDxtsA x coordinate of traffic signal upper left corner ytsDytsC y coordinate of traffic signal upper left corner xtsxtsAxtsBxtsCxtsDxtsA Vector of x coordinates of traffic signal corners ytsytsAytsBytsCytsDytsA Vector of y coordinates of traffic signal corners plotxtsyts k Plot outline of traffic signal hold on plotxtsgytsg o MarkerSize 40 MarkerEdgeColor none MarkerF aceColor g Visible visibleg plotxtsyytsy o MarkerSize 40 MarkerEdgeColor none MarkerF aceColor y Visible visibley plotxtsrytsr o MarkerSize 40 MarkerEdgeColor none MarkerF aceColor r Visible visibler axis square axis32 42 5 25 END PLACEMENT OF TRAFFIC SIGNALS Script file TrafficSignalCallm Script file that calls function TrafficSignalm to position traffic signal face and control which colors are lit TrafficSignal25 25 on on on 25 u 20 15 10 5 32 34 36 38 40 42 143 Function Workspaces Function mFiles create their own temporary workspaces for storage of variables created inside the function The variables are local to the function workspace meaning they are not accessible in the MATLAB Workspace or any other function workspace This explains why the same name can be used for a variable in a script file which is stored in the MATLAB Workspace and a variable created in a function mFile which is stored in the function39s workspace After the lnction executes the function workspace is deleted and the variable remains only in the MATLAB Workspace In the example below a scalar 39m39 and a lgtlt2 array 39pt are saved in the MATLAB Workspace and then passed to a function which uses the same variable names to create data for plotting a line Example 1431 function ortholinempt This function receives the slope m and a 1 by 2 array pt of x and y coordinates for a line It calculates the slope of the orthogonal line through the same point and then plots the original line in red and the one perpendicular to it in blue x0ptl Store 1st component of array pt in x0 y0pt2 Store 2nd component of array pt in yo xlinspacex05x052 Create 1 by 2 array x of end values yly0 mxx0 Find coresponding y values on line m1m Find slope of orthogonal line y2y0 mxx0 Find corresponding values on orthogonal line vx05x05y05y05 Set scale for plotting plotxyl r xy2 b Plot both lines axisv axis square The Command window statements are m2 pt10 40 ortholinempt In this example there are no variables created in the function 390 rtho l i ne and the function39s workspace would ordinarily be empty The numerical value of the MATLAB Workspace variable 39pt39 is read by the function but the array itself does not get copied into the function39s workspace The variable 39m39 on the other hand is copied into the function39s workspace because it is modi ed inside the function If necessary functions can share variables with the MATLAB Workspace if they are declared gl oba l in each workspace Example 1432 Script file cylcallm This file calls function cylm which returns the volume of a cylinder with radius R and height h The radius R is declared global in both files global R R4 Set value of global variable h10 Input passed to function cylm volsidecylh Call function cylm with input h Rhvolside function y1y2cylh This function is called by script file cylcallm and calculates the volume of a cylinder with radius R and height h The radius R is global in both files It also calculates the side of a cube with the same volume The volume is returned in output yl and the cube side is returned in output y2 0909090909 global R YlPihRA2 Cylinder volume YZY1A13 Side of cube Running the scn39pt le cylicall results in R 4 h 10 vol 5026548 side 79510 144 Function m Files and the MATLAB Search Path When MATLAB encounters name on the right hand side of an assignment in the MATLAB Command window or inside a script or function mFile it follows a hierarchial set of rules to determine what to do 1 It looks for name as a variable in the current workspace 2 It checls for a builtin mction name 3 It checks for a sub lnction name in the m File where n ame appeared 4 It checls if name m exists in the current directory 5 It checks if name m exists in any of the directories on the MATLAB search path The MATLAB search path is searched in the order in which it is speci ed When you create your own mFiles they should be stored in directories which have been added to the MATLAB search path Section 1 Introduction It is not uncommon to be working with two variables eg x and y that obey a fixed but unknown relationship y fx between them The function fx relating the two variables may not be as important as merely having the ability to generate by measurement or experimental observation numerical values x y where y fx For example the performance of a component may be related to a design parameter in a complicated fashion which is difficult to describe in mathematical terms In a controlled environment the lone parameter can be varied and the resulting performance monitored Suppose a consumer testing magazine is interested in the repair costs for a particular luxury passenger vehicle after its been in an accident The test procedure calls for headon collisions of the vehicle with a stationary object at various speeds The cost of restoring the vehicle to its precrash condition is determined Tabulated results after several collisions are given below Speed S Damages D mph 10 4500 20 19750 30 43500 40 55000 Table 11 Experimental Vehicle Crash Test Data The table represents a sample of four data points obtained from an unknown function D fS The experimental process is rather expensive and further experimentation might well be cost prohibitive The magazine would like to publish a graph for its readers to estimate the cost of repairing the same vehicle resulting from collisions occurring over a range of speeds One approach to approximating the unknown function D fS is illustrated in Figure l l The approximating function is the result of connecting the data points from Table 11 by straight line segments The dotted curve represents the actual unknown function D fS To estimate the damages from a 15 mph collision the approximating function is evaluated at this speed as shown The process of approximating a function x by another function and subsequent use to estimate numerical values of fx is referred to as interpolation The original function fx may be known in analytical form or possibly only through tabulated values that come from it The approximating function hereafter referred to as the interpolating functlon ls ordlnarlly mueh slmpler than the underlylng funeuon x 1n Flgure 1 1 the y x 6 X la 5 D S A 5 if 4 3 5 g 3 E 2 5 a 2 l 5 1 Approxrmatlng 4 12125 Fu etlon Estamated Repau n 5 Cost at 15 mph m u an 35 A 2 25 Speed mph Flgure 1 1 An Approxlm aung Funeuon For the Data m Table 1 1 erp atrng funetrdns are generally restneted to an lnterval eontarmng the data pornts used to generate them In 15 examp e estrmated repalr eosts should be ednfrne to speeds from 0 to 40 mph an addlnonal data pdrnt ean be assumed at S 0 mph D extrapolauon Cautron should be used when extrapolatrng slnce erroneous and even nonsenslcal results are posslble terpolatron of data representrng future predeuons ls not extrapolatron The two aphed m Flgure 1 2 represent hrstoneal data as well as forecasts of we economlc vanables Presumably sound economlc forecastlng methods were employed m the process of extrapolaung the future values Regardless The preeewrse 1mear rnterpolatrng funeuon ls a surtable apprommaung functlon g well behaved andth dt polnts Bylnlar e and e a a are reasonably spaeed then a eonneet the dots approach uslng preeewrse 1mear rnterpolatron wlll produce sausfaetory results 2000 WEED Social Security Receipts and Outlays 1600 MUD 39 39 g Oudays E 1200 V g mm 1 Reeerpts BBB BBB AUD 200 1a 5 20m 2005 201m 2015 2020 20 5 yeav Figure 1 2 Data Pomts from Funetaorrs Predretarrg Future Values As an example consider the ease where the function to be approxrrrrated is 1 t 1 a to 2 rad x x y m 7 sm x 0 1 00 0 8415 1 1 25 0 9490 2 1 50 0 9975 3 1 75 0 9840 4 2 00 0 9093 Table 1 2 Several Points From the Functionx sm x function sin r are shown in Figure 1 3 By observation rt appears that the approxrmatrrrg function is within afew pereerrt ofthe sine funetaor overthe errare interval P m Frgure 1 3 15 the estarhate a reasonably e1ose approxrrhataoh to the aetua1 value of m m2 oh the other hand 1f srh m2 was estarhaed rh t e hhe Jorhrhg porhts 1 0 8415 and 2 0 9093 then the resu1twou1o1 be they coordmate of pothzwh1h dAffErs srghr eahuy from the eorreet value 11 12 13 14 15 I113 17 1E 19 mad 2 Frgure13 Lmear 1nterpo1ataoh of x srh xto Estarhate srh m2 m graphmg sottware to p1ot known aha1ytaea1 fuhetaohs The fuhetaoh 15 eva1uateo1 at h e e Frgure 1 the porhts where the funcuon 15 eva1uateo1 e approxrrhatahg funcuon beeorhes smooth as beeorhe more dense The funcuon bemg approxrrhatedr e graphedr Preeewrse approxrrhatahg fuhetaohs are hot hrmted to the ease where a hhear funcuon approxrrhates the true funcuon between data porhts Later m the ehapter e eohsrder other po1yhorhra1 approxrrhatrohs to the uhder1yrhg funcuon x betwe eohseeutave dataporhts snaemn1 Spacmg u 5 1mm 1m an an 1 2 3 A 1 2 3 A snaemnu25 SpacngElEH 1mm 1m an an 1 2 3 4 1 2 3 A X X Flgure l 4 Graph of x ex Uslng Llnear Interpolau on Between Sampled Pornts 1n eontrast to the preeewrse approxlmanon of a funeuon rnterpolaung funeuons are frequently eornposed of allnear eornbrnatron of elementary funetrons o x x 01 1 Th M k 1x an gtnxar rx M x 1 1 gt x m h t F 1n tlus ease the rnterpolaung funetron ls a polynomlal denoted hymn where x2aruanaxazx2 my 12 en The nth orderTaylor Senes expmslon ofa functlon x mtxoducedln Chapter 1 ls Interpolatrng polynomlals are dlscussedln detarl m the remammg seeuons are several reasons why rnterpolatrng polynomlals are popular when rt There eornes to approxrrnauon of funetrons There are effluent algorrthrns for enaluatrng polynomlals when n ls large or the polynomlal rnust be eornputed nurnerous tunes wth different arguments Differentiation and integration of complex functions is often required In some instances the function is only available in tabular form Polynomial approximations of these functions are easily differentiated and integrated Exercises 1 The value of the Social Security trust fund at the end of several years is given in the table below Source US Social Security Administration 6 V O V D V Year Fund Value at End of Year in millions of dollars 1985 35842 1987 62149 1989 155063 1991 267849 1993 369322 1995 458502 Data for Problem 1 Plot the data points and the piecewise linear function connecting the data points on the same graph Estimate the fund39s value at the end of 1994 from a line connecting the following pairs of data points 1985 35842 and 1995 458502 1987 62149 and 1995 458502 1989 155063 and 1995 458502 1991 267849 and 1995 458502 1993 369322 and 1995458502 The fund39s value at the end of 1994 was 413460 million Calculate the true error for each estimate in Part b Use the piecewise linear interpolating function to estimate the fund39s balance at the end of June 1990 Use the last two data points to extrapolate the year and month that the trust fund exceeded 5 billion The actual occurrence was December 1996 US Census data at the end of selected decades is portrayed in bar chart form Source US Census Bureau a b Estimate the Us population in the year 1900 using a linear approximating function connecting data points for the years i 1880 and 1910 ii 1850 and 1940 iii 1820 and 1970 Us population in 1900 was 760 million Calculate the errors for the estimates obtained in Part b y 1 1 tt to estamate the years m whteh the U s popu1attoh reaehed 100 muhoh and 200 muhoh d Extrapo1ate the U s popu1attoh m 1980 and 1990 usmg the 1asttwo data pomts e Sketeh a smooth eurve through the rst four data pomts and use tt to mterpo1ate the US popu1auoh m the year 1840 Repeat the proeess for the 1ast four data pomts and estamate what the U s popu1attoh was m r e aetua1 popu1auohs were 17 1mt1hohm 1840 and 179 3 mt1hoh m 1960 2581 U s Pupmattun mtutuns 179m 1821 1851 man 191m 194m 197m yEav Graph for Problem 2 3 Solar radJauon is an tmportaut factor m eomputahg the eoo1mg1oadrequtremehts of e nthloner Measurements of so1ar radJauon through 18 m sheet glass facmg East were taker at dAfferent ttmes of the day at several 1atttuo1es or August 21 The E39 a Plot the data pomts wtth so1ar raduauoh oh the vemcal axts and lautude oh the honzohtal axts b P1ot the pteeewtse 1mear fuhetaoh for eaeh tame of day e Esttmate the so1ar radJatlon at 45 1atatuo1e at 9 am noon 3 pm and 6 pm d Plot the data points with solar radiation on the vertical axis and time of day on the horizontal axis e Plot the piecewise linear function for latitude f Estimate the solar radiation at 11 am for each latitude Latitude Standard Time 240 320 400 48D 56D 9 am 202 202 199 195 188 noon 71 70 68 65 61 3 pm 32 31 30 28 26 6 pm 9 9 10 11 11 Solar Radiation Heat Gain Through Glass Btuhr per ftz 4 Humidi ers are rated in terms of capacity ie gallons of water per day it can evaporate The moisture requirements of a tight house well insulated and an average house of different sizes are tabulated below The figures are based on a Winter day with outdoor drybulb temperature 200F 80 relative humidity an indoor drybulb temperature of 70 0F and 40 relative humidity Source Refrigeration and AirConditioning ARI Institute PrenticeHall VOLUME OF RESIDENCE TIGHT HOUSE AVERAGE HOUSE ft3 Gallons per Day Gallons per Day 8000 509 1017 16000 1018 2035 20000 1272 2544 28000 1781 3651 Table for Problem 38 Plot the data points for both houses on the same graph 2 V b Draw the piecewise linear approximating functions for estimating humidifier rating from the residence volume for both types of house c Estimate the size humidifier required for a tight house with 10000 ft3 of space d The humidifier in an average house must evaporate 315 gallons of water per day to maintain the design indoor conditions on a day when outdoor conditions are the design values Estimate the volume of the house e Are the answers to Part c and d exact values or simply estimates Explain 5 A digital sound le is obtained by sampling an analog audio signal using an analog to digital ND converter The digital le is played back by the reverse process using a digital to analog DA converter The following table lists the kb of memory required to store 1 second of audio for various sampling rates using two different ND converters It also shows the duration of audio that can be stored in 1 Mb of memory for different sampling rates 8 bit ND Converter 12 bit ND Converter Sampling Rate f kHz Length of Sound Stored in 1 Mb T sec Length of Sound Stored in 1 Mb T sec Memory Storage for 1 sec of Sound S kb Memory Storage for 1 sec of Sound S kb 4 4 20 40 250 8 50 40 25 80 125 25 125 20 40 O V Table for Problem 5 Plot the data points f s for both ND converters on the same graph Determine the equation of the interpolating polynomials for predicting memory storage requirements for 1 second of sound as a function of sampling rate for each type of converter Plot both polynomials on the same graph with the data points Plot the data points f t for both N D converters on a different graph Determine the equation of the interpolating polynomials for predicting the length of sound stored in 1 Mb as a function of sampling rate for each type of converter Plot both polynomials in Part e on the same graph with the data points of Part c A telephone conversation is to be digitally recorded Assume the sampling rate will be at twice the highest frequency component of the audio signal roughly 4 kHz Estimate the memory storage required to store a one minute conversation using an 8 bit ND converter An FM radio signal is to be digitally encoded using a 12 bit ND converter The frequency range is 20Hz l6kHz and the sampling is performed at 8 times the highest frequency component How many Mb of memory are needed to store a 3 minute song Is your answer to Part h the result of interpolation or extrapolation Is extrapolation reliable in this situation Explain Chapter 26 2D Graphics 261 The plot Function The basic function for producing 2D graphs is the MATLAB plot39 function It is generally called with two identically sized row or column vectors one for the set of abscissa values and the second for the ordinate axis The plotted points are ordered pairs of values one from the rst vector and one from the second vector For example to plot the points x1y1x2y2xn yn the arrays x c1x2 xn and y y1 y2 yn are created and then passed to the plot function for graphing Example 2611 x0 2 5 8 10 Create x vector y0 4 3 3 6 Create y vector P10txy Plot 00 24 53 83 106 When the plot39 command 393 issued within an mFile or from the Command window MATLAB opens a new Figure window draws a set of axes automatically scaled to t the data plots the points in the order in which they appear in the two vectors and then connects all the points by straight lines Numerical scales appear on both axes as do tick marks The blot39 command will clear the current Figure window if one is open before plotting the new data Several data sets can be plotted from one plot39 command by adding additional arguments Example 2612 xlinspace02pi500 usinx YCOSX zlsinxA2 z2cosxA2 zzl22 v0 2pi 11 11 plotxuxyxzlxz2xz axisv hold on ploto 2pi0 o39k39 Note if separate plot39 commands were issued ie plot x u plot x y etc only the last data set of points x1zlx2zzg002500would appear in the Figure window The arguments containing the data points to be plotted are not restricted to row vectors They may be rectangular arrays provided the number of rows and columns are the same The corresponding columns are paired and the data points plotted in the same way when the arguments are one dimensional row or column vect01s Example 2613 X18 916 1724 X is a 8 by 3 array Yrand83 Y is a 8 by 3 array of random numbers plotXY Plot 1st column of Y vs 1st column of X plot 2nd column of Y vs 2nd column of X Plot 3rd column of Y vs 3rd column of X x 1 9 17 2 10 18 3 11 19 4 12 20 5 13 21 6 14 22 7 15 23 8 16 24 y 07948 08757 02844 09568 07373 04692 05226 01365 00648 08801 00118 09883 01730 08939 05828 09797 01991 04235 02714 02987 05155 02523 06614 03340 1 I I I 09 08 07 06 05 04 03 02 01 0 I I I I 0 5 10 15 20 25 MATLAB also accepts a column vector as the rst argument and a recatangular array for the second argument as long as both have the same number of rows In this case each column of data in the array is paired with the column vector rst argument and there is one plot for each combination Example 2614 xlinspace01500 Create 500X1 column vector x YxxA2xA3xA4xA5 Create 500X5 array Y plotxY Plot columns of Y vs column vector x 39Plot39 can be called with a single argument say a column vector y of real numbers and the result is a graph of the elements of y vs the integers l 2 3 lengThW Example 2615 y05rand201 Create 20X1 column vector y P10ty Plot 1Y1 2Y2 20Y20 If y is complex the result is a plot of the imaginary parts of y vs its real parts Example 2615 wlogspace22500 sjW rho05 Dampling Ratio Wn10 Natural frequency Tf90sA2 2rhowns Wnwn Second Order System Transfer Function plotTf Polar Plot axis1 1 1 1 axis square hold on P10t0011 k391100 k xlabel Real Part of Tf ylabel Imaginary Part of Tf title Second Order System Polar PlotU Imaginary Part of Tf Second Order System Polar Plot 0 Real Part of Tf 262 Linestyles Markers and Colors Linestyles markers and colors of MATLAB plots are controllable by specifying a string argument immediately following the data The choices are Color Marker y yellow point m magenta 0 circle 3 cyan x xmark r red plus g green star b blue 5 square w white 1 diamond k black v triangle A triangle lt triangle gt triangle p pentagram h hexagram down up left right MATLAB uses a default color scheme to draw multiple curves on the same plot The order is easily determined as evidenced by the following example Example 2621 x01 Y1l 22 33 44 55 66 77 88 99 plotxY axis0 1 0 10 In the next example a set of data points is shown as blue diamond symbols and the points are connected using cubic splines drawn as a continuous red dotted curve Example 2621 x0021 Create x data yexpxrandl62 Create y data plotxy bd Plot data points as blue diamonds hold on xilinspace0l500 yiinterplxyxi spline plotxiyi r Plot cubic spline fit as red dotted curve The default line style is a solid line which will be used to connect the points with straight lines when there is no marker speci ed There is no default marker and none will be drawn unless one is speci ed However when a marker is speci ed for placement at the data points the markers are not connected by a straight line unless a line style is given 263 Plotting Styles Different background colors are available for the axes and the gure in which it is embedded The 39co l 0 rd e f 39 command controls this feature COLORDEF Set color defaults COLORDEF WHITE or COLORDEF BLACK changes the color defaults on the root so that subsequent figures produce plots with a white or black axes background color The figure background color is changed to be a shade of gray and many other defaults are changed so that there will be adequate contrast for most plots Note that when using a black background color on the axes the default color order for drawing line plots is different than when the default 39 colorde f white39 is in effect Example 2621 is run again with the 39 colordef black39 command to show the new ordering Example 2631 colordef black x0 Y11 22 33 44 55 66 77 88 99 plotxY axis0 1 0 10 title Colordef Black Plot Color Order 39Colordef Black39 Plot Color Order 264 Plot Grids Axes Box and Labels Horizontal and vertical grid lines can be superimposed on a graph by issuing a grid on command 39Box 0 f f 39 will remove the box which borders the axes Labels on the x and y axes and a title are obtained by enclosing text strings inside the 39xlabel 39 yl abel39 and Ltitle39 commands respectively Finally text can be placed at point xy on the plot when the 39textx y syntax is used The center left edge of the text is positioned at the point xy Example 2641 xlinspace010250 ylx5A2 5 y2y1 plotxylxy2 hold on x5 ylx5A2 5 y2y1 plotxy1 39 39 xy23939 box off grid on xlabel x ylabel y1 and y2 title Translated Parabola and its Reflection text447 Peak of yl text447 Min of y2 Translated Parabola and its Reflection y1 and y2 265 Customizing Plot Axes The axes on graphs can be customized in a number of ways The general procedure involves issuing a 39plot39 command followed by one or more 39axis39 commands The axis features under MATLAB s control are listed below AXIS Control axis scaling and appearance AXISXMIN XMAX YMIN YMAXl sets scaling for the x and yaxes on the current plot AXISXMIN XMAX YMIN YMAX ZMIN ZMAX sets the scaling for the x y nd zaxes on the current 3D lot V AXIS returns a row vector containing the scaling for the current plot If the current View is 2D V has four components if it is 3D V has six components AXIS AUTO returns the axis scaling to its default automatic mode where for each dimension limits are chosen based on the extents of all line surface patch and image children AXIS MANUAL freezes the scaling at the current limits so that if HOLD is turned on subsequent plots will use the same limits AXIS TIGHT sets the axis limits to the range of the data AXIS FILL sets the axis limits and PlotBoxAspectRatio so that the axis fills the position rectangle 39nice39 This option only has an effect if PlotBoxAspectRatioMode or DataAspectRatioMode are manual AXIS IJ puts MATLAB into its quotmatrixquot axes mode The coordinate system origin is at the upper left corner The i axis is vertical and is numbered from top to bottom The j axis is horizontal and is numbered from left to right AXIS XY puts MATLAB into its default quotCartesianquot axes mode coordinate system origin is at the lower left corner axis is horizontal and is numbered from left to right The y a to top tick mark in size This 39s is vertical and is numbered from bottom AXIS EQUAL sets the aspect ratio so that equal increments on the xy and zaxis are equal makes SPHERE25 look like a sphere instead of an ellipsoid AXIS IMAGE is the same as AXIS EQUAL except that the plot box fits tightly around the data AXIS SQUARE makes the current axis box square in size AXIS NORMAL restores the current axis box to full size and removes any restrictions on the scaling of the units This undoes the effects of AXIS SQUARE and AXIS EQUAL AXIS VIS3D freezes aspect ratio properties to enable rotation of 3D objects and overrides stretchtofill AXIS OEE turns off all axis labeling tick marks and background AXIS ON turns axis labeling tick marks and background back on Examples 2651 through 2655 illustrate the use of some axes control features In the first example a circle and a square are plotted without any subsequent axis control commands Neither shape appears correct because of the screen s 43 aspect ratio The scales are also set automatically to contain the plotted data points Example 2651 h300 k500 r200 Set circle parameters xilinspace100500250 yuik sqrtrA2 xihA2 Find upper halfcircle y values ylik sqrtrA2 xihA2 Find lower halfcircle y values plotxiyui b xiyli b Plot circle hold on xa200xb200xc400xd400 ya400yb600yc600yd400 xxa xb xc xd xa Fix x locations of square s corners yya yb yc yd ya Fix y locations of square s corners plotxy Plot square 300 I I I I I 100 150 200 250 300 350 400 450 500 In the second example the MATLAB code is identical to Example 2651 except for the addition of the last command 39axi 5 equal39 which forces the aspect ratio to 11 changing the appearance of the the circle and square h300 k500 r200 Set circle parameters xilinspace100500250 yuik sqrtrA2 xihA2 Find upper halfcircle y values ylik sqrtrA2 xihA2 Find lower halfcircle y values plotxiyui b xiyli b Plot circle hold on xa200xb200xc400xd400 ya400yb600yc600yd400 xxa xb xc xd xa Fix x locations of square s corners yya yb yc yd ya Fix y locations of square s corners plotxy Plot square axis equal I I I I I I I 50 100 150 200 250 300 350 400 450 500 550 Note that the plot is still enclosed in a rectangle and the scaling is still chosen by default to assure the entire circle is Visible In the next example the scales along both axes are set with the appropriate axis scaling command The gure is an equilateral triangle which appears to be an isosceles triangle Example 2653 L2 Equilateral triangle with sides of length 2L xA3 xBxAL xCxA2L Set x corners of equilateral triangle yAl yByAsqrt3L yCyA Set y corners of equilateral triangle v0 10 0 5 vxmin xmax ymin ymax plotxy Plot equilateral triangle hold on axisv Set x and y scales to values in vector v xdata0 xA xBxC2 ydata0 yA yByC2 ppolyfitxdata ydata 2 Fit quadratic thru xdata ydata xi inspace0xBxC2100 yipolyval pxi plotxdata ydata k xi yi r Making the border of the plot square using axi 5 square39 does not correct the problem as illustrated in the next example Example 2654 Equilateral triangle with sides of length 2L xB xAL xCxA2L Set x corners of equilateral triangle axisv Set x and y scales to values in vector v xdata0 xA xBxC2 Ydata0 YA YByC2 ppolyfitxdata ydata 2 Fit quadratic thru xdata ydata xilinspace0xBxC2100 yipolyval pxi plotxdata ydata k xi yi r axis square To render the equilateral triangle with equal sides on the screen the axis square command is replaced by axis equal as shown in Example 2655 Example 2655 L2 Equilateral triangle with sides of length 2L xA3 xBxAL xCxA2L Set x corners of equilateral triangle yAl yByAsqrt3L yCyA Set y corners of equilateral triangle axisv Set x and y scales to values in vector v xdata0 xA xBxC2 ydata0 yA yByC2 ppolyfitxdata ydata 2 Fit quadratic thru xdata ydata xilinspace0xBxC2100 yipolyval pxi plotxdata ydata k xi yi r axis equal Note the ymin0 and ymax5 values are no longer in effect 266 Multiple Plots As we have already seen its possible to include multiple plots on the same set of axes The hold on39 command retains previous plots in the current set of axes and remains in effect until a 39hold of f 39 command is issued Subsequent plots are drawn on an empty set of axes and cleared prior to the next plot until hold on 39 is reissued The current hold state ie on or o ff can be ascertained using the 39ishold 39 command which retums a 1 if hold on39 is in effect and 0 ifnot Example 2661 ishold if ishold0 Check if hold is off hold on ishold end if ishold for i1z5 xlinspace0lrand500 yx i plotxy Color rand rand rand end for i1z5 hold off ishold g 5 m H Ol O 267 Multiple Figures More than one Figure window can be used for plotting sets of data The 39figure39 command issued in the Command window or from within an mFile creates a new Figure window with an integer handle to identify it The last Figure window created is the active or current gure To make a different Figure window active either click on it with the mouse or issue the command 39figure h 39 where h is the handle of the Figure window to be activated The current Figure window is closed by a mouse click in the upper right hand comer or as a result of the MATLAB command cl 0 s e 39 Any existing Figure window can be closed using close h 39 where h is the handle of the Figure window to be closed All Figure windows are closed with close all 39 Finally 39cl f39 will erase the contents of the current Figure window without deleting the Figure window The following example contains a script mFile which opens three Figure windows and plots different sets of data in each It should be run from the Command window and not the m Book Example 2671 Script file uniformaccm This file will accept initial conditions for a uniform acceleration trajectory and plot in separate Figure Windows the dispacement velocity and acceleration profiles 09090909 close all Close all figures a0input Enter constant acceleration in ftsec sq Accept uniform acceleration value x0input Enter initial displacement in ft Accept initial displacement v0input Enter initial velocity in ftsec Accept initial v city xa02 v0 x0 Set displacement polynomial coefficient vector vpolyderx Find velocity polynomial coefficient vector apolyderv Find acceleration polynomial vector tilinspace010100 xipolyvalxti vipolyvalvti aipolyvalati plottixi Plot displacement data points in Figure 1 xlabel t sec ylabel xt ft title xt vs tU figure Open Figure 2 plottivi Plot velocity data points in Figure 2 xlabel t sec ylabel vt ftsec title vt vs tU figure Open Figure 3 plottiai Plot acceleration data points in Figure 3 xlabel t sec ylabel at ftsec sq title at vs t After running the script le Figure 3 is on the screen and Figures 1 and 2 are minimized on the task bar Clicking on Figure l in the task bar will make it the active Figure window The plot of displacement Xt vs t in Figure l is shown below for values of a02 X0100 and v05 Xt vs t l l l l l l l l l Uniform acceleration a 4 ftsec sq 200 t Y 150 100 I I I I I I I I I 0 1 2 3 4 5 6 7 8 9 10 268 Subplots More than one set of axes can be postioned in a Figure window In fact the window can be segmented so that several rows and columns of axes will appear simultaneously To con gure the current Figure window for an mxn arrangement of subplots the MATLAB command is 39s ubplot m n p 39 where p refers to the active subplot ie the area where the next data set is to be plotted The numbering of subplots is from le to right starting in the top row The following example simulates 20 outcomes of an experiment It graphs all outcomes in the first plot upper left comer and then continues to graph subsets of the outcomes in succeeding plots Example 2681 x10randl20 subplot22l Activate upper left corner of 2X2 subplot Window plotltx3939 axis1 20 0 10 xlabel n Trial Number ylabel 0utcomes All subplot222 Activate upper right corner of 2X2 subplot Window xlt5indicesfindxlt5 xlt5xxlt5indices plotxlt5indicesxlt5 U axis1 20 0 10 xlabel n Trial Number ylabel 0utcomes lt 5 subplot223 Activate lower left corner of 2X2 subplot Window xgt5indicesfindxgt5 xgt5xxgt5indices plotxgt5indicesxgt5 U axis1 20 0 10 xlabel n Trial Number ylabel 0utcomes gt 5 subplot224 Activate lower right corner of 2X2 subplot Window xbtwnindicesfindxgt25 amp xlt75 xbtwnxxbtwnindices plotxbtwnindicesxbtwn axis1 20 0 10 xlabel n Trial Number ylabel 0utcomes gt 25 amp lt75W Outcomes All Outcomes gt 5 n 439 1 qu 4quot 5 10 15 20 n Trial Number l arr 5 10 15 20 Trial Number Outcomes lt 5 Outcomes gt 25 amp lt75 Iquot y 15 20 Trial Number n 10 15 20 Trial Number 269 Interactive Plotting Tools The 39legend command places string labels in a legend box to help distinguish between the different curves within a set of axes Example 2691 includes plots of several trigonometric functions and a legend to help identify each one Example 2691 xlinspace02pi500 ysinsinx ycoscosx ytantanx V0 2pi 3 3 plotxysinxycosxytan axisv legend sinx cosx tanx hold on plot0 2pi0 0 k Plot x axis xlabel x rad 3 l l I l l I sinx cosx 2 tanx The legend box can be dragged to another location if the default position is unsatisfactory A zoom39 function in MATLAB allows you to zoom in or out inside a Figure window It can be activated and deactivated by typing 39zoom on and zoom off39 in the Command window or the statements can be included in an mFile When 39zoom is on left and right mouse clicks zoom in and out from the mouse point by a factor of 2 The 39zoom39 feature can also be turned on by selecting the magnifying glass icons in the gure window tool bar Other features of zooming include the capability of zooming in either the x or y directions only as well as specifying a zoom factor Type help zoom39 to learn more about zooming Example 2692 x100011 ysinxx subplot2ll plot xy xlabel x ylabel sinxx subplot2l2 Plot XIY axis1 1 0 2 xlabel x ylabel sinxx zoom3 Zoom in by 2 and 2 again Warning Divide by zero 25 E I 03 1 00 00 04 02 0 02 04 00 00 1 X 12 11 25 1 E I 09 0398 I I I I I I I I 02 015 01 005 0 005 01 015 02 In the above example the zoom3 command zooms in on both axes by a factor of 2 two times That is the xaxis goes from 11 to 0505 to 025025 and the y axis changes from 02 to 05l5 to 075l25 Also note the 39hole39 at x0 in the magni ed plot and the NaN waming MATLAB has amdamentary dlgmzlng capablllty whlnh ls desmbed below 5mm Graphlca npu from mouse X31 GINPUTINI gens N P lnts from ch courdlnaces n length n veeeere x and y e eur xx 7 GINPLVT geehers er unllmlced number k sed uncll ehe reeurn ey ds pres quot w othe d l d Flgure wlndow The glnput39 funedon ls used to extract the coordinates of points along Gemlm Blvd to the North A spllne function ls tted to the points and men plotted as an overlay to the lm e Example 36 10 3 meadt cmusaexlal xmmw ls Read Image as was me lxl Display my Cnnecl lxm cu rdmau of ram 31an mad any road in srnne evaluallnn 1 uemua n2357 muss H4579 n59uu H5544 anus EIE7EE 1EI277 11755 13135 14EI25 14524 15425 15515 1 uemua n 5953 H595 EI6EI4 H6615 n77ss H534D n5225 755u H735 H752 n925u LEIEES 115m 12995 2611 Specialized 2D Plots Data can be plotted on logarithmic scales by using either semilogx39 39semilogy39 or loglog39 instead ofthe conventional plot39 command Example l6lll xlinspace0510 y310Ax2 semilogyxy b xy r39 The function used to generate the data points was y 3 10 2 Starting with the plotted data points the function can be estimated using the end points as follows yAlOBx thru x0 y3 andx5y950 log3logAB0 3 A3 log950log3B5 3 B Wwojoo1 The 39comet function provides an eyeappealing way of plotting connecting data points The following example must be run from the Command window or a script le to demonstrate the come t39 plot feature Example 26112 Script file CometPlotm This file plots two comet style graphs xlinspace02pi500 Ysinx cometxy Plot sine function pause2 close1 thetalinspace02pi500 r2sin2theta xypol2cartthetar cometxy Plot rose curve pause2 close1 thetalinspace010pi500 r2theta xypol2cartthetar cometxy Plot spiral A stacked area plot permits a cumulative function and the contribution of each component to be graphed with the areas between the components filled in For example fx xfzxf3x x2xl f2x xi j x sin727x the next example plots the stacked components of fx using the MATLAB area39 Jnction Example 26113 xlinspace02pi500 f13x1 f2xA2 f320sinxA2 Ff1 f2 f3 areaxF xlabel x ylabel fxflxf2xf3x legend flx3xl f2xxA2 f3x20sinxsinx 2 f2xx2 f3x205inxsinx f1xf2xf3x fa A polygon described by a set of vertices can be drawn and lled in using the f i ll 39 mction An example follows Example 26114 xlsqrt22 y10 x5sqrt2 y5sqrt22 mly5ylx5x1 y2sqrt26 x2xly2ylml y4sqrt23 x4xly4y1ml x3x5y3y2 x6x2y6y5 x7x1 y72sqrt23 x8x7x6x7 y8y6 x90 y9y8 x10x1x4x1 Y10y4 x110y11y3 x12x8y12y11 xx1 x2 x3 x4 x5 x6 x7 x8 x9 x10 x11 X12 X11 yy1 y2 Y3 Y4 Y5 y6 Y7 YB Y9 y10 yll y12 yl axis equal fillxy r When two different sets of data require different yscales the solution is to use MATL s 39plotyy39 function Example 26115 xlinspace 0 2pi250 Ysinx z5expx plotyyxyxz xlabel x ylabel y legend sinx 5expx Since 39plotyy39 actually creates two sets of axes inside a Figure window its ill adVised to issue any commands which alter the appearance of axes since they will affect only one of the two sets of axes and produce unintended results In fact in the preVious example the legend command was not recognized for the x2 axes Try adding the command 39axis tight alter the plotyy39 command to the code in Example 26ll5 and observe the results 4000 2000 There are several other specialized 2D plots available in MATLAB Bar charts stem plots and polar plots are shown in the following examples Example 26116 xl 12 Monthlytemp6010randl 12 30sinpi x 16 bar x Monthlytemp Plot bar chart xlabel Month Number ylabel Ave Monthly Temperature 80 7O 60 50 4O 30 Ave Monthly Temperamre Example 26117 xl 12 Monthlytemp6010randl 12 30sinpi x16 bar3 x Monthlytemp Plot 3D bar chart xlabel Month Number ylabel Ave Monthly Temperature Ave Montth Temperature Month Number Example 26118 k0z20 ek4 06 quotk 5 08 quotk stemkek Plot discrete data xlabel Discrete Time k ylabel Discrete Error Discrete Error 039 46 TTTIW 8 1o 12 14 Discrete Time k Example 26119 thetalinspace0 2pi 250 r12sin theta polarthetar title Polar Plot of a Limacon Polar Plot of a Limacon 90 3 Section 3 Newton DividedDifference Interpolating Polynomials The standard form for representing an nth order interpolating polynomial is straightforward There are however other representations of nth order polynomials which on the surface may seem a bit more unwieldy but require less manipulation to arrive at One such form is the subject of this section Given a set of n1 data points xi yi i 0 l 2 n where the x are all different and the yi are sampled from an underlying function y x the nth order interpolating polynomial can be expressed as 1395 b0 bloc 960bzx xox xibnx xox xi39quotquotx xm 31 Before we consider how to determine the coef cients observe that Equation 31 is in fact an nth order polynomial as evidenced by the last term which includes the highest power of x namely x All powers of x are present in Equation 31 despite the fact that the coef cients of say x0 x1 x2etc are not as obvious as when the polynomial is expressed in stande form see Equation 21 Since there are n1 independent coef cients b available 139 0 l 2 n we can be assured there is at most an nth order polynomial that passes through the given data points The rationale for selecting the analytical form in Equation 31 will be apparent after we look at a few simple examples A Given x0 yo x1 y1 where yo lm and y1 xl The rst order polynomial for the case when n 1 through the two data points is f1x b0 b1x x0 32 Substitution of the two data points x0 yo x1 y1 into Equation 32 gives f1x0b0 b1x0 x0 33 99bob1xl xo 33a Solving for b0 and b1 yields b0 f1xO 34 bl 34a x1 x0 By design the interpolating function f1x and the actual function fx from which the data points were obtained are equal at x x0 and x x1 see Figure 22 As a result b0 and b1 can be expressed in terms of the given data be fx0 35 b1fx1fxo x1 x0 Example 3 1 The monthly payment on a 30yr mortgage of 100000 for two different annual interest rates is given in Table 31 below Use an interpolation formula in the form of Equation 32 to estimate the monthly payment corresponding to an interest rate of 825 per year Data Point Number Annual Interest Rate Monthly Payment k ik Ak z k 0 7 66530 1 10 87757 Table 31 Monthly Payments for 100000 30yr Mortgage with Different Interest Rates Two Data Points The first order interpolating polynomial is written f1l39bob1i io 36 where b0 fl39o f7 66530 7076 b fz391 fz390 f10 f7 87757 6653 1 1391 1390 10 7 10 7 The estimated monthly payment is therefore f1 825 6653 7076825 7 75368 B GiVen 95070 bel xzy2 Where yo xo yl xl and yz xz The second order polynomial is written f2xb0b1xx0b2xxoxx1 where b0 b1 and b are determined using the same procedure employed in the previous example f2x0b0b1xox0b2xoxoxox1 38 f2x1b0b1x1x0b2x1xox1x1 38a f2x2 b0 b1x2 x0b2x2 xox2 x1 38b Solving for b0 b1 and b2 gives b0 fxo 3399 bl fx1fxo x1 x0 fx fxfxmm x x x x b 2 1 1 39b x2 x0 Once again each fzxi was replaced by fxi so the coefficients could be expressed in terms of known quantities There are de nite advantages to representing interpolating polynomials in the nonstandard form In the previous example there is a sharp contrast in the method of solving for the three coefficients between the standard form representation of a second order polynomial in Equation 27 and the nonstandard form of Equation 37 In standard form a0 a1 and a are obtained by solving the system of simultaneous equations given in Equation 210 On the other hand b0 b1 and b are obtained sequentially as indicated in Equations 39 39a and 39b Given the choice a sequential solution is usually preferable The savings in computation is real It is faster to evaluate explicit formulas for a set of coefficients one at a time than it is to implement a matrixbased solution where the coefficients are obtained in a parallel fashion ie solution to a system of simultaneous equations This argument may be less convincing when using calculators or computers with programs designed to solve simultaneous equations Nonetheless the solutions are obtained in fundamentally different ways The second advantage of the nonstandard representation is more compelling The expressions for b0 and b1 in Equations 39 and 39a are the same as in Equation 35 Why is this important Suppose you just completed the process of finding an interpolating polynomial for a given set of data points There may be some doubt in your mind concerning the accuracy of results based on the use of this polynomial Later on there will be a discussion of quantitative methods for approximating the errors inherent in interpolation More data points may be required to reduce the estimated errors Incorporating additional data points is easier with the nonstandard representation for the interpolating polynomials This is because each additional data point requires the computation of a single coefficient for the new term in the polynomial This is illustrated in the following example which extends the results obtained in Example 31 Example 3 2 Suppose we obtain two additional data points in the previous example dealing with the estimation of mortgage payments The new data points correspond to 8 and 9 loans Use one of the two additional points to obtain a second order interpolating polynomial Estimate the monthly payment for an 825 loan We must choose one of the two new data points to find the second order interpolating polynomial fzi Our intuition suggests the data point for an 8 loan is the wiser choice because its closer to the point where the estimate is required ie 825 The second order polynomial is obtained from the data points corresponding to k 0 l and 2 in the following table It is given in Equation 310 Data point Number Annual Interest Rate Monthly Payment k ik Ak ik 0 7 66530 1 10 87757 2 8 73376 3 9 80462 Table 32 Monthly Payments for 100000 30yr Mortgage with Different Interest Rates 4 Data Points f2ib0b1ii0b2iioii1 Coefficients b0 and b1 were determined in Example 31 using the first two data points The remaining coefficient b2 is obtained from Equation 3 as fan filfil fio 12 11 11 10 I l 2 0 73376 87757 87757 66530 8 10 10 7 7 1148 The polynomial f21 is therefore f2 1 6653 70761 1011481 101 11 311 and the estimate of payments for an 825 mortgage is now f2 825 6653 7076825 7 1148825 7825 10 75124 The previous two examples are somewhat academic in nature Quite obviously the loan officer at the bank will not be interested in your calculations to estimate the monthly payments He or she has access to the functionA f1 from which you obtained several data points in the Sunday paper Let s see how close your two estimates were to the correct answer Figure 31 includes graphs of the true function f1 as well as the interpolating polynomials f11 and f21 The upper plot includes the data points as well It is clear from looking at the upper plot that either interpolating function f11 or fz1 will provide accurate estimates of the true monthly payment function f1 over the entire range from 1 7 to 1 10 and then some The lower plot is an enlargement of each graph in the region about 1 825 the interest rate under consideration Here we see that the quadratic interpolating polynomial fz1 and the real function f1 are virtually indistinguishable You may be interested in knowing that the true function A f1 is given by M2001200quot A P 312 fa 1 X200quot 1 where P is the mortgage amount 71 is the loan period in months and 1 is the annual percent interest rate Perhaps you are surprised at how close the graphs of the low order interpolating polynomials f11 and fz1 are to the true function f1 which based on its analytical form in Equation 312 does not appear to resemble a polynomial function 5 Many Paymem A 5789111111213141 Annua1 1n1cycs1 Raw 1 37m lt 2751 Emu S fun E73m EDS 81 E15 82 325 83 Annua1 1n1cycs1 Raw 1 835 EA 845 as Frgurc 3 1 Fxrst and Sccoud Order Interpolanng Polynomxals x x and the Truc Functionx Table 3 3 xllustxates mc accuracy ouc can expect when usmg low order polyuorrua1 mterpolauon over the range ofmterest rates under cousrdcrauou x ltgt I x AC ET 13920 ET 7 5 699 21 700 68 7147 699 25 70 04 8 5 768 91 77143 V2 52 768 86 0 05 9 5 840 85 84219 7134 840 76 0 09 Table 3 3 Errors m Fxrst and Sccoud Order Interpolauon for Examples 3 1 and 3 2 We can have an addmonal data pomt that cou1o1 be used to mcrease the order of the mterpolatmg po1yuorrua1 The formulas for hrghcr o der coef ments bccoruc uwrc1o1 Fortunately there 15 a framework for deter mmmg the coef ments 5 x 0 1 n f a general mh order mterpolatmg po1yrrorrua1 as gven m Equauor 3 1 The general approach 15 based on the cva1uauor of fxmte dwrdcd dAfferences Gwen a set of n1 data pomts rx x 0 1 2 n the mte dwrdcd dAfferences of various orders are defined in Equations 313 313e Note how the diVided differences are obtained recursively from two diVided differences of order one less 0thorder fxl fx 313 1storder fxlx W i j 313a xl x fog mp 3m xl x 2ndorder fxlxxk W i j k 313c xzxk fx mp fx x 2 M 313d xzxr nthorder fxnxn1mxpx0 fxnxn1x2x1 fxnilxn2x1x0 3136 xn xo When there are 3 data points 71 2 the diVided differences fx0 fx1 x0 and fxz x1 x0 are identical to the coefficients b0 b1 and b2 see Equations 314 314h The coefficient b0 of the interpolating polynomial fzx is numerically equal to the first of the three zero order finite diVided differences ie the one that depends on the data point 9 fx0 The coefficient b1 is equal to the first of the two first order finite diVided differences ie the one that depends on the data points 9 fx0 and x1 fx1 Finally the coefficient b2 is equal to the first and only second order finite diVided difference ie the one that requires all three data points x0fx0 x1 fx1 and x2 xz In general with n 1 data points there are 7 1 zero order diVided differences 7 first order diVided differences nl second order diVided differences etc up to one nth order diVided difference The first computed diVided difference of order quot139quot is equal to the coefficient b 139 0 1 2 n flxol fx0 314 b0 314a xlaxo M 3140 x1 x0 fx1fxo 31 1196259615960 l 172 fx2fx1 x1 x0 flxzax1fx1x0 x2 x0 H x2 x1 x 2 0 fx1fxo x1 x0 314d 3l4e 3143 l 314g 314h Table 34 is a helpful aid in remembering how to compute the finite divided differences The coefficients bi 139 0 l 2 n of the interpolating polynomial in Equation 31 are the first finite divideddifference entries in their respective columns The polynomial is referred to as the Newton divideddifference interpolating polynomial The arrows indicate which finite diVided differences are used to calculate the higher order ones between its first and last argument Always remember the denominator of any diVided difference is the difference alll fn fll l f2 l 960 flxol x1 x1ltfx1 x0lgt xz x1 x0 fx2 x1 x2 fxzlt fx3 x2 x1 flx3 x2 x393 flJ39C3l x3 flx4 93963 x2 xnZ fxnz fxn139 xn2 xnla xn 2 x2 x0 xquot xn1quot39x2 x1 fxn xn1 xnZ x1 xo xn1 fxn1 fxn xn1xnz gt xquot xnl xn flxn Table 34 Table of Finite DiVided Differences for a Function x Known at Discrete Data Points Example 3 3 Torquespeed data for an electric motor is given in the rst two columns of the table below Find the equation of the Newton diVideddifference interpolating polynomial that passes through each data point and use it to estimate the torque at 1800 rpm Speed 0 Torque T1 f1 f2 f3 f4 rpm x 1000 ft1b 05 31 b0 6 b1 10 28 2 b2 8 6667 b3 15 24 12 6 b4 20 5333 2 0 14 4 24 25 2 Table 35 Data Points and Finite DiVided Differences for Example 36 The notation f1 fz f3 and represent nite diVided differences of order 1 through 4 respectively The nite diVided differences and the coef cients b0 b1 b2 b3 and b4 of the interpolating polynomial are shown in the table The result is Ta 0 b1a w0b2w a0a w1b3w w0w w1w wz b4a 60060 60160 a2a a3 315 Ta 31 6 a 05 2 a 05a 1 6667 a 05a 1a 15 6 w 05a 1a 15a 2 316 T18 31 6 18 05 2 18 0518 1 6667 18 0518 118 15 6 18 0518 118 1518 2 T18 187 The true function T f0 is not available to assess the accuracy of this estimate Conceivably one could derive a structure or analytical form of the function based on established engineering principles and scienti c laws A number of physical parameters ie constants would likewise have to be known before the function could be used for evaluation of motor performance Rarely is this ever attempted in situations where all that is required is a reliable estimate of how a particular device component or system will perform Conversely a model of the torque motor based on scientific principles is far more useful to an engineer designing a motor to satisfy specific design criteria There are situations where a combination of empirical and scientific modeling are used in conjunction with each other Imagine a situation where a theoretical model is known however it may be of such complexity that data points are expensive to obtain In this situation a minimum number of data points can be obtained from solution of equations differential and algebraic comprising the scientific model The data points are then used as a basis to generate an empirical model for interpolation Little has been said up to this point about the ordering of the data points primarily because an nth order polynomial passing through n1 data points is unique and the ordering is irrelevant Changing the order of the data points will affect the representation of the polynomial however the polynomial itself has not changed The selection of which n1 data points to draw from a larger sample to obtain an nth order polynomial is a different matter In this case some thought should be given as to which n1 points produce the best interpolating polynomial Example 34 illustrates these points Example 34 Consider the data in Table 32 Find the third order Newton divideddifference interpolating polynomial f3i containing the 4 data points using different orderings of the data points Estimate the function value f825 The data from Table 32 and the divided differences are given in Table 36 k ik Akfik f1 f2 f3 0 7 66530 b0 6846 b1 1 8 73376 1200 1 7086 005167 b3 2 9 80462 1045 7295 3 10 87757 Table 36 Finite Divided Differences and Newton Interpolating Polynomial Coefficients for Data in Example 32 The coefficients b0 b1 b2 and b3 are read directly from the table The resulting third order Newton divided difference interpolating polynomial is given in Equation 318and the polynomial is then used for the interpolation of f825 f3z39 b0b1i iob2i ioi ilb3i ioi ili iz 317 f3 139 653 6846z39 7 12139 7z39 8 005 167139 7z39 8z39 9 318 f3 825 653 6846825 7 12825 7825 8 005167825 7825 8825 9 75126 In Table 37 the data points are reordered and the process of finding the interpolating polynomial is repeated k ik Akfik f1 f2 f3 0 7 66530 b0 70757 b1 1 10 87757 1148 b2 71905 00515b3 2 8 73376 1045 70860 3 9 80462 Table 37 Finite Divided Differences and Newton Interpolating Polynomial Coefficients for Data in Example 32 with Reordering of Data Points The new polynomial based on reordered data points and the interpolated value are f3 139 653 70757z39 7 1148139 7z39 10 00515139 7z39 10z39 8 319 f3825 653 70757825 7 1148825 7825 10 00515825 7825 10825 8 75126 As expected the estimated values are the same in both cases It39s a simple matter to show that the polynomial expressions in Equations 318 and 319 are in fact the same third order polynomial Restricting the interpolating polynomial to be 2quotd order will produce as many different 2quotd order polynomials as there are ways to select groups of three points from a total of four data points In general for a speci c value of 139 referred to as the interpolant each interpolating polynomial fzi will yield a different estimate of t One such estimate of f825 has already been computed using the first three points in Table 32 The optimal choice of data points for determining the second order interpolating polynomial should be the set of points closest to the interpolant This will be proven later Generally speaking when Newton divideddifference polynomials x are used for interpolation to estimate values of a function fx we should expect a difference between x and the true function value fx In the case of linear interpolation the error incurred R1x is the difference between x and f1x and it varies over the interval x0 S x S x1 as shown in Figure 32 By definition R1x satisfies f x 1200 R106 320 Solving for R1x R1xfxf1x 321 We should not expect to implement Equation 321 when fx is unknown If there was some way to estimate fx we could obtain an approximation for the error term or remainder R1x as its sometimes referred to From Figure 32 observe that fzx could be used for that purpose That is R106 f2 x 1200 322 The right hand side of Equation 322 is nothing more than the second order term of the polynomial fzx This is easily verified by looking at Equations 32 and 37 Thus Rm z bloc x0x x 323 However in order to use Equation 323 an additional data point xz xz is required to obtain b2 How good is the approximation of R1x from Equation 322 or equivalently Equation 323 Comparison of Equations 321 and 322 reveals that the approximation of the error term R1x depends on how close the function x and the interpolating polynomial fzx are Figure 32 illustrates this point graphically If the estuhate of Rm would be exactfor any value ofx A humeheal example wlll help clarlfy the process of approxlmaung the ehth A Kn X2 xl thure 3 2 Estlmauhg the errorR1x lh Llheat Interpolation Example 3 5 The solutaoh to the dlfferenual equataoh ytytl wheh y0 0 ts yt le equot12 0 Let T be the tame lt takes for the solutaoh yt to reach the vain A e where 0 g A g l l e y Several polhts T A on the solutaoh yt are tabulated below A low orderpolynomlal ls neededfonnterpolanon ofT for aglven value of A 03567 03 16094 08 05108 04 23026 09 Table 38 Data Points for Finding Interpolating Polynomial to Estimate T Givenl The underlined data points were selected to determine a third order Newton diVideddifference polynomial which can be used for interpolation over the interval The diVided difference table is shown below z39 A TifAi f1 f2 f3 0 00 00000 11890 1 03 03567 13275 21183 47745 2 07 12040 76676 3 09 23026 40238 4 05 06931 Table 39 DiVided Difference Table for Finding f3A and A From the table the third order interpolating polynomial f3A is f3 A 1189A 13275AA 03 47745AA 03A 07 324 Suppose we wish to estimate f045 the time required for the solution to reach 045 using f3045 The result is f3045 118904513275045045 03 47745045045 03045 07 05441 Choosing the additional data point 05 06931 because of its proximity to the interpolant value 045 allows us to nd the 43911 order interpolating polynomial A Table 39 includes the additional finite diVided differences calculated from the new data point The result is f A 1189A 13275AA 03 47745AA 03A 07 76675AA 03A 07A 09 325 14 and the improved estimate of f045 is therefore f4 045 f3045 76675AA 03A 07A 09 05441 76675045045 03045 07045 09 05441 00582 06023 Based on reasoning analogous to that used to obtain R1x in Equation 322 the error term R3045 is approximated as R3045 2 f4 045 f3045 2 00582 It is easily shown that the true function relating T andA is given by T fA ln1 A 326 Figure 33 shows two graphs Each one contains the complete set of data points from Table 38 and the true function fA The top graph contains the third order polynomial f3A and the four data points used to find it are shown with an asterisk The lower graph contains the fourth order polynomial A The same four data points used to find f3A and the one additional point used to find A are shown as asterisks The true error R3045 is also shown in Figure 33 As expected R3A is zero when A is any of the four data point values The exact value for R3045 is obtained as the difference between f045 andf3045 The result is R3045 f045 f3045 1n1 045 05441 00537 which should be compared to the estimated value of 00582 previously computed Keep in mind that f3A A and R3A are all sensitive to the choice of data points selected from Table 38 Knowing a priori that f045 is to be estimated would dictate a different choice of data points as the basis for determining an interpolating functlon to appronrrnate A The 4 elosest pornts to tlne rnterpolant value 045 are 0 30 3567 0 40 5108 0 50 6931 and 0 60 9163 25 2 Flgure 3 3 Thlrd and Fourth order Interpolanng Polynomlals and tlne True Funen on for Datarn Table 3 8 s pornt rt ls neeessary to lntxoduce a new At th functlon x r rm slmllar to a seeond order dvrded dlfference de ned as fmm xlmrlrlxrxnl 327 n 3 27a However unllke tlne seeond order nrte dlvlded dlfferences prevlously eneountered rts rst argument x ls treated as a v able Of eourse onee x assumes a n al In A x1x es an or nary seeond order dvrded dlfference The rrnportanee of x n n ls rts relauonshlp to Mr whlch ls presented here and let as an exerclse for later R106 x x0 x x1fxx1xol 328 If fx is a slow varying function over the interval x0 x1 then fx x1 x0 will likewise vary slowly over the same interval Under these conditions for any value of x in the interval x0 x1 fxx1xol fxzxpxol 329 Replacing fx x1 x0 with its approximation fxz x1 x0 in Equation 328 gives R106 x xox x1fx2xlxol 330 which is identical to Equation 323 It was previously stated that if fx happens to be a quadratic polynomial then the estimate of R1x in Equation 322 would be exact What does this imply about the function fx x1 x0 under the same conditions We can now generalize the preceding analysis about the error in Newton divided difference interpolating polynomials With n1 data points x fx 139 0 12 n from a function x the Newton divideddifference interpolating polynomial x that passes through each point is related to the function x in the following way fxfnxRnx 331 The error term Rnx accounts for the difference between the true function x and the interpolated value x See Equation 320 for the case when 711 Equation 328 is a special case 711 of a general expression for Rnx which follows Rnx xxoxx1xxz39quot39x xn1x xfxxnxn1x1xol 332 The nlst order finite divideddifference term f x Equation 332 cannot be evaluated at the interpolant value of x unless the function fx is known Instead an estimate of f x x1x0 is possible if an additional data x1x0 in xwxm xwx point xn1fxn1 is available The estimate is fx xrl lelquotquotquotxl x0 z fxn1 xn xn15 quot5x15x0 The right hand side of Equation 333 is of course bn1 the coefficient of the nlst order term of 1x Using this approximation in Equation 332 provides a means of estimating the error term Rnx The result is Rnxaxx0xx1xx2nnx xn1x xnfxn15xn5xrl1quotquotquotx1 x0 Rn x bn1x xox x1x x2 39 39 39 39x mex 96 335 17 The estimate of Rnx is therefore the highest order term of the nlst degree divideddifference polynomial 1x ie the new term added to x made possible by the acquisition of an additional data point In the case of linear interpolation 111 the estimate of the error term R1x from Equation 335 reduces to the expression given in Equation 323 Equation 332 for Rnx is important for several reasons Imagine a situation where exhibits little or no variation over the interval corresponding to the n1 data points xifx 1 0 l 2 n In the extreme case where fx x1x0 is constant Equations 333 and 334 are no longer approximations and the true value of fx is computed from Equation 331 In this case the interpolating polynomial fn1x and x are identical The following example illustrates this situation xnxn1 Example 3 6 Three data points xifxi 1 0 l 2 from the function fx x3 l are given in the table below A second order divideddifference interpolating polynomial fzx is required The finite divided differences are also given in Table 310 i xi xi flll f21 0 1 0 13 1 3 26 10 63 2 6 215 Table 310 Divideddifference table for Data in Example 36 The second order divideddifference interpolating polynomial is f2 x 13x l10x lx 3 336 If we use Equation 336 for interpolation the error term Rzx can be approximated after we obtain an additional data point from fx Instead three additional data points x fx 1 3 4 5 are calculated from the cubic function x and the complete listing of finite divided differences is given in Table 38 We will show that the estimate of the error term Rzx is exact regardless of which of the three additional points is used to estimate Rzx Moreover we will see that the third order divided difference fx x2 x1 x0 is constant Choosing x3 fx3 5124 as the new data point the third order divided difference fx3 x2 x1 x0 1 From Equation 334 Rzx is approximated by 18 R200 5 x39x0x39x1x39x2fx3x2x1x0 R2 x z xlx3x6gtl 338 i x fXi f1 f2 f3 d fs 0 1 0 13 1 3 26 10 63 1 2 6 215 14 0 91 1 0 3 5 124 13 0 39 1 4 2 7 11 28 5 4 63 Table 3 ll Divided Difference Table 6 Data Points From A Cubic Function fx x3l Using f2x for interpolation at x 35 f2 35 1335 l1035 l35 3 45 and estimating the error Rz35 from Equation 338 R2 35 a 35 l35 335 6 r 1 t 3125 How good is this estimate of R235 We can answer this question because the function fx is known The true error is R235 foofies 353 11 45 3125 19 The estimate of Rz35 is exact and f35 will be computed exactly from f335 An expected result since f3x and x are different representations of the same third order polynomial function In this example the estimate of Rzx is exact regardless of where the interpolation occurs The reason is simple We learned that Rnx in Equation 334 is exact whenever the nlst divided difference f x x1x0 is constant ie independent of x This is precisely what happened in this example namely fx x2 x1 x0 is constant See if you can show why this is true whenever x is a cubic polynomial xwxm Table 311 includes three additional data points despite the fact that only one extra point is needed to approximate Rzx Observe that all three third order divided differences fx3 x2 x1 x0 fx4 x3 x2 x1 and x x4 x3 x2 are equal Let s reconsider this example from a different perspective Suppose the function x is in unknown and we have obtained the six data points x xi 139 0 l 2 3 4 5 from it given in Table 311 Unaware of the true function we set out to compute all the finite divided differences intending to employ a Newton divideddifference interpolating polynomial of order as high as 5 if necessary In the process of completing the table however we notice that the column of third order divided differences has the same numerical value The unknown function x is exhibiting a property characteristic of third order polynomials namely fx x2 x1 x0 is constant However we need to be careful before concluding that the function x is a cubic polynomial It is correct to say there is a third order polynomial that contains our six data points It is given by f3x13x 1 10x lx 3 1x lx 3x 6 339 Assuming that f3x and the true function x are one in the same is not justified simply because f3x passes through the six data points sampled from x The higher order interpolating polynomials x and f5x reduce to f3x because b4 and b5 are both zero If f3x and x are the same function then subsequent points obtained from x would produce the same result as we saw in Table 311 namely all third order finite divided differences would equal one If there exists a function x not a cubic polynomial that happens to pass through four five or even all six data points of the function x3 1 then an additional data point from x would result in a different value of the third order divided difference Can you think of such a function Figure 34 is a graph of two possible functions x along with f3x Note that one of the possible x functions passes through four points of f3x and the other passes through five points of f3x The first function was generated in MATLAB by specifying 20 4 pornts from g and one addluonal pornt not on f3x The polyflt funeuon of order4 x poly lt the seeond funeuon x whlch passes through 5 pornts of g and one addlnonal pornt not on f3x x through 5 pornts of g x through 4 polnts of g Flgure 3 4 Two Noanublc Functlons x wrth Pornts Common to g In summary mferences eoneernrng the strueture of x eannot be made solely on The rmportanee of the prevlous dlscusslon ls not whether a functlon x ls aetually a polynomlal Rarel ls the underlymg funeuon a true polynomlal What s r p trs th le dl de tYerenees reveals m a tan way lf polynomlal wlll be surtable for rnterpola o If a p r ar eolumn of dlvlde dl erene s ls s stant en an rnterpolaung pol omlal emsts that a g 8 aeeurately matehes the glven data pornts Furthermore tlus polynomlal than the nu order polynomlal whlch ls guaranteed to pass through the n1 data pornts as the followmg example rllustrates Example 3 7 Measurements of coffee temperature taken at 2 minute intervals are given in the table below Find the lowest order polynomial suitable for interpolation of the function relating temperature and elapsed time Time I Temp Ti f1 f2 f3 f4 f5 fst min 0F 0 0 212 225 1 2 167 1875 150 0 1250 2 4 137 1125 00130 105 00208 00021 3 6 116 1000 00078 000028 65 00833 00013 4 8 103 0500 00052 45 00417 5 10 94 0250 35 6 12 87 Table 312 DividedDifference Table For Data Points in Example 37 Results of the divided differences in Table 312 suggest that we should consider a third order polynomial for interpolation To be sure lets graph the third order polynomial based on the data points ti Ti 139 0 1 2 3 This polynomial is f3t 212 225t l875tt 2 0125tt 2t 4 340 and is shown in Figure 35 as ff I Not exactly what we expected Using f3t in Equation 340 for interpolation would result in significant errors for t in the interval from 8 to 12 minutes Of course if you like your coffee at 120 F the interpolating polynomial will do just fine There are several options we can explore at this point A higher order interpolating polynomial is certainly one of them The 4m 5m and 6th order polynomials can be obtained and then graphed to assure that no anomalies or pronounced uctuations are present as a result of using higher order polynomials 22 Ifwe are rntent on usmg a thrrd po1ynomra1 forrnterpo1atron we ean combme the four thrrd order fxmte dvrded ddfferences t 1211 11314114 ta tts 1413 12 and ma ts to 13 and use the aver e value m p1aee of tg 1 tr tn m formrngsa The resu1trng thrrd order rnterpo1atmg po1ynomr shown m Frgure 35 as 5 15 an 0 a1 rmprovement over the one grven m Equanon 3 4 f ZerZZSI18750727129070074 ffc2127225t1879072 11115771073074 WEIEI Data Pomts for 3 tn 0 Tan 212 ED 70 167 m2 137 EU 70 11s 2 A E E 1D 1 mm Frgure 3 5 Thrrd order 1nterpo1anng Po1ynomra1 From F1rst4 Data Pornts m Table 3 12 ere 15 yet another alternauv which m nts consideration The decision to e e dd dedddff nesmT o by a etor apprommate1 70 0208 to 701250 1mplymg the thrrd order ddfferencefuncuon A m m anve1y eonstant as mmally assumed Ihrs on W111 be sensrnve to the ehoree of base 112 tr ay not be re1 means the po1ynomra1 f3t used for mterpo1atr pornts usedto de ne rt e better t for f3t ean be obtarned Ordenng othe data pornts 15 rrre1evant only when the entrre set 15 usedto determrne the umque po1ynomra1whreh passes througq eaeh pornt If we rntendto use O for rnterpo1anon an here between 0 and 1 enset 1 e1ude the rst and1ast measured data pornt m the determrnatron of f3t Table 39T M V 2 mmutes then 1tmakes 23 Rational approximations were used in the MATLAB execution to minimize roundoff errors in the calculations The new interpolating polynomial is f3t212 16ttt 6 tt 6t 12 341 Time I Temp Ti f1 f2 f3 f4 f5 fst min 0F 0 0 212 16 l 6 116 6772 1296 ll8 2 12 87 1724 52304 254 796 74608 3 4 137 916 l92 1346080 l72 596 l768 4 8 103 1312 0 323 596 5 2 167 3748 738 6 10 94 Table 313 DividedDifference Table for Example 37 with Reordered Data Points and the graph in Figure 36 is more to our liking Notice the smaller variation in the column of third order divideddifferences in Table 313 compared to Table 312 There are other choices for the 4 base points that determine f3t but we shouldn39t expect much if any improvement over the interpolating polynomial in Equation 341 After all the true function Tt for the cooling coffee is certainly not a cubic polynomial We now explore the relationship between the nite divided differences and the derivatives of the interpolating polynomials The order of a polynomial function is implied by the highest order derivative which is a constant For example the first derivative of a first order polynomial is constant the second derivative of a quadratic is constant etc In other words the nth derivative of an nth order polynomial is constant and the n1 st derivative is zero In Equations 342 and 343 below fx is an nth order polynomial 24 MK viir eorm lt3 42 air and dm M a dxm fltxgt 0 3 43 31 212e1s621zesexxzzesze 12 1am 1am 2141 12m mm Data Pomts for 3ng 10 0 73900 212 z 12 73902 87 EU 134 739tg137 13 1mm 1n Frgure 3 s Thrrd Order Interpolatmg Po1yrrorma1 A er Reordenng 0mm Pomts Supposethe nth orderpolynomxal x passes through me n1 pomts x 1a x 012 n We have already seen that the mh order dwxdedrdAfference funeuor 15 constant See Example 3 s and Table 311 That rs n 1W constant 3 44 Subsmuung Mr from Equauon 3 32 mto Equation 3 31 yre1ds xIXX39Xnx39xlx39xz X Xr 1X Xvavx1v vxxvxn 3 45 The mterpolatmg polynomxal x W111 be 1denuca1 to the nth order polynomxal x andtherefore the terrorism m Equation 3 45 mustbe zero Thatxs xxoxx1xx239quotx xn1x xnfxxmxn1axpxol 0 346 Since Equation 346 must hold for all x we conclude that the nlst divided difference function satis es fxxnxn1x1x00 347 In summary for an nth order polynomial function fx f 1 x constant 342 fxxn1 x1 x0 constant 344 f 116 0 343 fxxnxn1x1x0 0 347 The constants in Equations 342 and 344 are not equal however the fact that both fquotx and fx xn1 x1 x0 are constant is the important point Equations 342 through 344 along with Equation 347 suggest the existence of a relationship between finite divided differences and derivatives when fx is an nth order polynomial Exploring this further consider the case where n 2 and the polynomial function fx is given by fx a0 alx a2x2 348 Its left as an exercise to show that the second divideddifference function fx x1 x0 is constant as is f2x What about the first divideddifference function fx x0 and the first derivative f x Figure 37 shows a quadratic function x and a graphical interpretation of the first divided difference function fx x0 keeping in mind that x is variable For a given x fx x0 will be the slope of the line segment connecting x0 and x There is a point C located between x0 and x where the first derivative f C and fx x0 are equal Its fairly straightforward to show that C is given by x0x 2 349 What happens if x in Figure 37 is no longer a quadratic function As long as fx is continuous and differentiable over x0 x the Mean Value Theorem from Differential Calculus guarantees there will be at least one C where the first derivative f C and fx x0 are equal 26 wnen x ls no longer an nth order polynornral none ofthe denvanve funenons are rdenneally al zero for x However e Lh order dlvl edrddffere new x E andthemh denvanve fo are stlllrelatedln aslmllarw relanonslnp ls presented here Wthout proof It applles for any value of 1 rm ampmfg where exx xn n mu x 2x1au Lar 2z 5 Slope x rn l 51092 R0 4D 2D El 1 2 3 A 5 E 7 E 9 1 z r Flgure 37 FlrstDlvldeerlfferenceFunctlon AME andLheFlrstDemanve f kx wny ls thls useful7 Recall that when a eolurnn of dlvlded dlfferences m a dlvldedrddfference table ls relatlvely constant we expect a polyn conespondlng order to be well sulted for lnterpolanon Refer to Table 3 13 Where we eoneluded that a thlrd order lnterpolanng polynomlal should work well beeause all four rd 4 a a a H mm thrd denvatlve of the underlymg functlon T0 m Example 3 7 must also be relatlvely constant In fact we e approxlmate it uslng the average 3quot order dlvlded dlfference on me le slde oquuanon 3 50 Dolng thls ornlal of the Wm 3 12vtlvtn13v12vtltov13v121f51513 351 27 e functh ay n The 3 50 67 to Incidentally the third derivative of the cubic interpolating polynomial in Equation 341 the one based on reordering of the data points to include the entire cooling period of 12 minutes is easily computed to be l3 Equation 350 guarantees that when the computed nth order nite divided differences are small the nth derivative of the underlying function is likewise small over the entire interval associated with the data points Hence the order of the interpolating polynomial is established with con dence When more data points are available than needed for a given polynomial our intuition suggested use of the data points closest to the interpolant See Examples 32 and 35 Our intuition has a solid analytic foundation The error term Rnx resulting from the use of an nth order interpolating polynomial was given in Equation 332 and repeated below Rnx xxoxx1xxzquotquotx xn1x xnflxaxnaanauaxpxol 332 The n1 st divided difference f x function fx is unknown Our only alternative is to choose the n1 base points x fx 139 012n from the larger set of data points to minimize the product of the linear factors in the expression above x1x0 is unobtainable when the xmx In our cooling coffee problem Example 37 if we require an estimate of T3 the cooled temperature after 3 minutes the interpolating polynomial in Equation 340 based on the initial four data points in Table 311 is the logical one to use In this case the product 3 to 3 t1 3 t2 3 I3 is minimized and hopefully the term R33 given below as well R33 3t03t13t23t3 3at39t29t19t0 352 One look at Figure 35 reminds us that the same polynomial is not the one to use for estimating T10 We have seen how its possible to estimate the error Rnx incurred when estimating x with an nth order interpolating polynomial fnx All that was needed was an additional point and the use of Equation 334 Unfortunately the application of Equation 334 can underestimate Rnx by a significant amount 28 What would be nice to have is an upper bound on Rnx Let s see if this is possible From Equation 332 and Equation 350 with 7 replaced by n1 Rnxxx0xx1x xnw where g e xxnxn1x0 353 n1 Suppose M is an upper bound not necessarily the maximum on the absolute value of the nlst derivative fquot1x over the interval containing the data points It follows from Equation 353 that an upper bound on the absolute value of the error is M I w me4 Ja m 354 Equation 354 is of little use when fx is available only in tabular form as sampled data points Nonetheless there may be times when an upper bound for the nlst derivative fn1x is obtainable Of course if the upper bound M happens to be the maximum absolute value of the nlst derivative fquot1x over the interval the upper bound on l Rnx l in Equation 354 is even better Example 38 A uniform cable hanging under its own weight is suspended between two poles 500 ft apart as shown in Figure 38 The cable sags 100 ft at its lowest point Measurements of the cable height were taken at several points and included in Table 314 xm ym xm ym 250 4280 50 3318 200 3909 100 3434 100 3434 200 3909 50 3318 250 4280 Table 314 Measured Points Along Suspended Cable Due to symmetry only the point x 0 y 328 and the last four points were considered These ve points were reordered the divided differences computed and entered in Table 3 15 29 100a x x y m 2 f f4 0 0 328 0 1 250 428 0 04810 60000ee007 2 50 3318 0 001 4 2 0000 e009 0 3940 8 0000 e007 3 200 390 9 0 00162 0 47 100 43 4 3 Table 3 15 Dmdedertferenee Table For Example 3 8 lung at me magnxtudes of me hrgher order dwxded dAerrences r1 appears that can be approxxmated very well by a second order polynomxal 2 gt4 Loo are true funcuony From Table 31129 15 Mr32804r0001o2rre250 3 55 30 A graph of the quadratre rhterpolatrhg polynomlal x based on the rst three porhts m Table 3 l5 and the rrght half of e eable pro lex are shown m Flgure 3 9 The polynomlal x passes very elose to the two addluonal pomts mm and row as peeted ex 45D Klryl ADD f2x32804x000162xx7250 gt Interpolauhg Polynomlal Based on m Karyn My mm W e L rezzaeush m eatehary my yo myu Kzryz Flgure 3 9 Cable Pro lex and Quadraue Interpolatrhg Polynomlal x An upper bound or the error R2x when usmg quadraue res some knowledge of the behavlor of the thud dehvauvefJ x w the thud dehvauve from Equauoh 3 50 usmg n x2 n m m Table 3 15 oh the elt ah slde ofthe equauoh to glve rhterpolatroh requl e ear approxuhate lo 3 x3xzxan for 09950 3 56 3 6gtlt10 7 36MBquot Suppose we wlsh to estuhate the eable helght 125 lt from rts lowest porht r e 125 Uslng Equatroh 3 55 glves f2125 352 69 ft From Equauoh 3 53 wth h 31 R2025 357 125 x0 125 x1125 x2 ag Approximating f6g by the value obtained in Equation 356 gives R2025 125 0125 250125 5036 x 10quot 422 Ordinarily there is no way of computing an upper bound on the error Rnx since we have no way of determining M an upper bound on fquot1x in Equation 354 However in this example the cable pro le function fx shown in Figure 39 assumes a form known as a catenary given by fx 328 cosh 358 Knowing the true function we can compute M the maximum value of the third derivative over the interval corresponding to the first three data points shown below f 3 96 An upper bound on the Equation 354 is R2125 R203491 The work is 32813 cosh i 359 dx 328 l x s1nh 360 3282 328 lMax f3x for 0 s x s 250 361 L 3282 250 sinh 328 77908 x 10 6 interpolation error Rzl25 using M 77908 x 10396 in s 125 x0 125 x1l25 x2 M 362 s 125 0125 250125 50 77908x10 s 913 32 We can compare this upper bound on the interpolation error Rzl25 with the true error because the function fx is known Doing so ET fx f2x 363 f125f2125 Lcosh 2 35269 328 328 058 The magnitude of the true error lETl 058 is of course less then the estimate of 422 obtained by approximating the 3rd derivative f3x using the third order nite divided difference fx3xz x1 x0 in Equation 356 It is also less than the upper bound 913 computed in Equation 362 which was based on the maximum 3rd derivative of the catenary function over the entire interval 0 S x S 250 33 l N Exercises A tutoring service has kept records of performance on a standardized test and the number of days students attend their review classes The performance rating Y represents the per cent improvement in the test score students attain after taking the exam a second time X is the number of attendance days in the review class x0 x1 x2 x3 x4 X Attendance Days 1 25 5 65 9 Y Improvement yo yl yz y3 y4 2 5 l l 14 17 Data Points for Problem 1 Plot the data points 2 V b Prepare a complete table of nite divided differences 0 Estimate 4 the improvement in one s score after 4 days of attending review classes assuming y x is the true function relating X and Y Base your result on a third order Newton divideddifference interpolating polynomial f3x d Plot f3x on the same graph with the data points e Estimate the error in the interpolated value f4 Use the extra data point ie the one not used in Part c to obtain your answer f Repeat Part e with x 0 y 0 as the additional data point Compare results of Parts e and f g Suppose the last point in the table was x4 0 y4 0 Repeat Parts a though e with x 3 instead of 4 The following table contains the x y coordinates of several points along the catenary shown in Figure 38 The variable 3 in the last column is the length of cable from the lowest point x 0 y 328 to the point x y i xi y Si 0 0 3280 0 1 50 3318 502 2 100 3434 1016 3 150 3629 1553 4 200 3909 2126 5 250 4280 2749 Data Points for Problems 2 and 3 34 E 4 UN VV C9 3 00 VV Plot the points x 3 139 0 l 2 3 4 5 Find the Newton first order interpolating polynomial f1x through the end points x0 30 and x5 35 Plot the linear function on the same graph with the data points Use the interpolating function f1x to estimate the length of cable required to span a horizontal distance of 150 ft starting from the lowest point on the cable 0 328 Suppose the true function relating 3 and x is 3 x Estimate the true error R1150 f150 f1150 by using the additional data point 200 2126 Find the second order Newton interpolating polynomial fzx based on the two end points x0 30 and x5 35 and the additional data point x4 34 Plot it on the same graph with f1x and the data points The true function x is 3 csinh Plot it on the same graph with f1x fzx C and the data points Calculate the true error R1150 Find an upper bound on the true error R1150 Refer to the same table of data points used in Problem 2 UN VV 0 g h Plot the points y 3 139 0 l 2 3 4 5 Find f1y the Newton first order interpolating polynomial through the end points yo 30 and y5 35 Plot the linear function on the same graph with the data points Use the interpolating function f1y to estimate the length of cable required to span a vertical distance of 349 ft sta1ting at the lowest point on the cable 0 328 Suppose the true function relating 3 and y is 3 y Estimate the true error R13629 3629 f13629 by using the additional data point 3909 2126 Find the second order Newton interpolating polynomial fzx based on the two end points yo 30 and y5 35 and the additional data point y4 34 Plot it on the same graph with f1y and the data points The true function y is 3 wlyz 3282 Plot it on the same graph with f1y f2y and the data points Calculate the true error R13629 Find an upper bound on the true error R13629 The region from sea level to approximately 36000 ft is called the troposphere where the temperature of the standard atmosphere varies linearly with altitude 35 a Describe the true variation T Th if the temperature at sea level h 0 ft is b O V Q V D V C9 51869 0R and 44743 0R at an altitude of 20000 ft Air pressure and the speed of sound vary with altitude in the troposphere The following table contains several data points re ecting the variation of these quantities with altitude Altitude h ft Pressure p psi Speed of Sound v ft sec 0 147 11164 5000 122 10971 10000 101 10774 20000 676 10369 30000 437 99485 36000 331 96875 Atmospheric Properties for Problems 4 and 5 Source Fundamentals of Flight R Shevell PrenticeHall 1989 An interpolating polynomial is required for estimating p over the entire range of altitudes Reorder the data points if necessary calculate the complete set of nite divided differences and enter the values in a finite divided difference table Determine the lowest order Newton divideddifference polynomial fnU l suitable for interpolation of p Do not choose the full fifth order polynomial ie n lt 5 Plot the data points hi p 139 012345 along with the interpolating polynomial p h where 1 S n S 4 on the same graph Use fn25000 to estimate the atmospheric pressure in the troposphere at an altitude of 25000 ft ie f25000 where p h is the true relationship Select one of the data points from the table that was not used to determine fnU l and use it to estimate the error in 25000 The true function relating p in psi and h in ft is p f0 p0 R where g gravitational constant at sea level 3217 ftsecz R gas constant for air 1718 ftlbslug OR a atmospheric temperature gradient in the troposphere ie slope of linear function T Th p0 atmospheric pressure at sea level psi 36 05 T0 atmospheric temperature at sea level OR Plot the function p h for 0 S h S 36000 on the same graph with the data points hivpi 139 012345 and the interpolating polynomialp h l S n S 4 g Calculate the true error f25000 25000 Refer to the same table of data points used in Problem 4 a An interpolating polynomials is required for estimating v over the entire range of altitudes Reorder the data points if necessary calculate the complete set of nite divided differences and enter the values in a finite divided difference table Determine the lowest order Newton divideddifference polynomial fnU l suitable for interpolation of v Do not choose the full fifth order polynomial ie n lt 5 b Plot the data points h v 139 012345 along with the interpolating polynomial v h where l S n S 4 on the same graph c Estimate the speed of sound in the troposphere at an altitude of 15000 ft ie fl5000 where v h is the true relationship d Select one of the data points from the table that was not used to determine fnU l and use it to estimate the error in 15000 e The true speed of sound at 15000 ft altitude is 10574 ftsec Calculate the true errorf15000 15000 A tank is filled with water to a certain level H0 A valve is opened and the time T required for the tank to empty completely is recorded The following data was obtained Initial Height Emptying Time H0 ft Tmin 50 589 37 507 28 441 10 264 5 187 Data for Problem 6 a Plot the data points with T as the dependent variable b Find the third order Newton divideddifference interpolating polynomial f3H0 suitable for estimating the emptying time when the tank is initially filled with somewhere between 20 and 30 ft of water Plot f3H0 on the same graph 37 gt1 O V Estimate the time it takes for the tank to empty when there is initially 25 ft of water in it Q V Estimate the error in your answer to Part c by using the remaining data point D V Estimate the error in your answer to Part c by using an additional data point of 0 ft 0 min Compare the answers to Parts d and e and comment on the results C9 The actual function for calculating the height of uid in the tank at any time is Z t 2A Ht Hf c for OStS H0 2A 0 where H 0 is the initial height of water in the tank ft A is the crosssectional area of the tank ft2 0 is a constant speci c to the tank ft3minft12 Plot the graph of T Ho for a tank withA 100 ft2 and c 12 ft3mimft1z g Calculate an upper bound on the absolute value of R325 25 f325 the true error in your answer to Part c 7 V Calculate the true error R325 The dynamics of physical systems are often modeled by second order differential equations of the form d d WWI Kw gym abt Kahlua where ut is the input and yt is the response When the input ut l for t 2 0 the output yt is called the unit step response and is given by l ytK 1 Ze 39sin ll Zant9 1 where 9cos391 and 0lt ltl The parameters of the system on C and K are called the natural frequency damping ratio and steadystate gain respectively A typical response is shown when on l radsec K l and 0 lt C lt l The percent overshoot PO refers to the amount by which the peak value M P of the time response overshoots the final value of K The P0 for various values of C are tabulated below a Graph the response yt with C one of the values in the table and verify the percent overshoot of the response is in agreement with the tabulated PO value 0quot lradsec andK l 38 b P1ot the data porhts gP o x 01 e Fmd the seeohd order Newton ddvxdedrddfference rhterpolatrhg polyhormal fag through the data porhts comespondmg to z 0 0 3 and 0 9 Plot the polyhorma1 on the same ga h as the data porhts and eorhrheht on how well you beheve rt approrrrhates the true fuhetroh P o er d Use f2 to estrrhateo 2 and0 8 e Estrrhate the true errors R20 2 u 2 exam 2 and Raw 8 U 8 exam 8 by orht ooPo 948 0 Calculate the true errors R20 2 and Raw 8 g Use the addmonal data porht m Part e to nd a thrrd order rhterpolatrhg polyhorma1 f3 P1otrt on the same graph wrth the data porhts and2O 4w h The true fuhetror P o g rs grveh by P o 1002W Plot 1 on the same graph wrth the data porhts audthe two rhterpolatrhg polyhorrua1s Frgure for Problem 3 7 0 01 02 03 04 053906 07 08 09 Po 100730527372254163948460152015 39 Table for Problem 37 8 Two points on a circle of 1 mile radius are shown in the diagram below The initial point 0 is xed and the second point P is located at coordinates xy The distance from O to P along the circle is D1 and the straight line distance is D2 The difference A D1 D2 is of interest Values of D1 and D were measured for various points along the circle The results are summarized in tabular form D1 POW O A Figure for Problem 38 1 xi miles yi miles D miles Dy miles A miles 0 000 00000 00000 00000 00000 1 025 06614 07227 07071 00156 2 050 08660 10472 10000 00472 3 075 09682 13181 12247 00934 4 100 10000 15708 14142 01566 5 125 09682 18235 15811 02423 6 150 08660 20944 17321 03623 7 175 06614 24189 18708 05480 8 200 00000 31416 20000 11416 Data Points for Problem 38 s V Plot the data points xi 4 1 0 1 8 U V Find the lowest order Newton interpolating polynomial that ts the data points reasonably well Plot it on the same graph Plot the full 8Lh order Newton interpolating polynomial f8x on the same graph 33 Plot the data points 11 4 1 0 1 4 on a new graph 40 gt0 D V Is it possible to find a Newton interpolating polynomial from the data points 01 A 1 0 l 4 to adequately approximate the true function A fyy Explain why the data points used to find the interpolating polynomial are restricted to a subset ofthose in the table ie only y corresponding to 0 S x S l Ifyou can find an acceptable interpolating polynomial plot it on the same graph Ifnot explain why CD Show that the true function A fxx is given by A cos 1l x 2x12 Plot the true function A fxx on the graph used for Parts a b and c Do either of the interpolating polynomials in Part b or c approximate the true function reasonably well 00 V h Find the equation of the true function A fyy Use it to generate several more points 01 4 Repeat Parts d and e using the expanded data set I VV Discuss what steps are necessary to find an interpolating polynomial to approximate A fyy for 0 Sy S l and l S x S 2 A digital filter is used to extract information from an analog signal corrupted with noise The analog signal is sampled to generate discrete data which is then numerically r J by a quot 39 39 r or a r 39 quot A digital signal processor DSP optimized to perform the required calculations The sampling rate is determined by the frequencies of the noise and signal components of the unfiltered wave form A common type of digital filter is called an infinite impulse response HR filter The order of an IIR filter is related to its noise rejection capabilities The amount of numerical computations increases with the order of the filter A DSP chip with U0 ND conversion and additional memory was tested to determine the maximum sampling rate achievable to implement IIR filters of various orders The results are summarized in a table 1 IIR Filter Order Maximum Sampling Rate n fmax kHz 0 l 1700 l 12 5 50 2 24 270 3 36 200 4 48 150 5 60 130 Table for Problem 39 a Plot the data points 11 fmax 1 015 4l Find the 5th order Newton interpolating polynomial f5n through the data points 6 V c In a particular application where the signal and noise frequency characteristics require a sampling rate of 400 kHz is it possible to use a 153911 order IIR lter using the tested DSP hardware Base your answer on the interpolating polynomial found in Part b d Use f5n to estimate the maximum sampling rate for a 301h order IIR digital lter Find the 2quotd order Newton interpolating polynomial f2n through the data points corresponding to lter orders of 12 24 and 36 Repeat Part d using fzn Compare fz30 and f530 and discuss which one is likely to be closer to the true value Two additional data points become available ie 20 320 and 50140 Which one would you choose to estimate the true error R530 30 f530 where fmax n is the true function relating the lter order n and the maximum sampling rate fmax Explain D V C9 g Choose the appropriate new data point and calculate the estimate of R530 A signalized intersection is the subject of traf c delay study Records taken over the same time period for a duration of several weeks indicated the major ow averages 1200 vehicles per hour from both directions and the minor ow averages 600 vehicles per hour entering the intersection from both directions The stopped time delay D refers to the time a vehicle is stopped on its way through the intersection Measurements were taken of the stopped time delay of every car through the intersection during the time period of interest Vehicles entering the intersection when the signal is green are assigned a zero delay since they do not stop The average delay per vehicle Dave is then calculated The cycle time C of the signal was xed at 90 seconds The portion of cycle time when the major ow receives a green light is designated GC The traf c engineer is trying to determine the optimum GC split during the time period of interest The following data resulted from the study 239 GC Dave secvehicle 0 055 326 1 060 322 2 050 346 3 065 330 4 045 408 5 070 355 42 Table of Data Points for Problem 10 s V Plot the data points GC Dave 139 01 5 U V Calculate the complete set of nite divided differences and enter them in a table Comment on the suitability of using Newton interpolating polynomials GC n 2 3 4 5 to approximate the underlying function Dave fG C Plot the Newton interpolating polynomials GC for n 2 3 4 and 5 es Estimate the optimum cycle GC ratio and corresponding minimum average delay for each polynomial in Part c Given two data points 9 x0 and x1fx1 from a function x demonstrate the validity of Equation 328 which states R1x xx0xx1fxax1xol 43 Section 5 Splines Up to tlus pornt we have restneted the appronrrnatrng funeuons that were used for rnterpolauon to a slngle polynomlal wth u conespondlng to the set of data pornts The results ean b number of reasons related to ether the order of the polynomlal or the general pro le of the data pornts When the order ls too low or too hrgh relauve to the vanablllty and el Flgure 5 1 shows a set of runs data polynomlal f2x through the mlddl andtwo polynomlal fed through all the data pornts suggests x approaehes a lrrnrtrng value m pornts from a functlon x a seeond order end data pornts and nally the elghth order A vlsual rnspeeuon of the data pornts an asyrnptoue fashlon ngh order Data Pornts For Deterrnrnrng f2x 01481951108 9988 Flgure 5 1 Example of a Funeu on x Not Well Surted for Polynomlal Approxlm an on Another approaeh to obtarn an appronrrnaung functlon for rnterpolauon ls based ornrnaung functlon n a e ewlse ashlon Infantw lgure 1 1 m Seeuon 1 rllustrate e s of preeewrse llnear segrnents to de ne the appronrrnaung funeuon Thls ls equivalent to using multiple first order interpolating polynomials restricted to intervals between adjacent data points Suppose we wish to approximate a function x de ned in tabular form by a set of data points x fx 139 012n with lower order polynomials each valid for a limited portion of the interval x0 S x S xquot These lower order polynomials are referred to as spline functions The term originates from drafting where a smooth curve is drawn through a set of control points using a exible strip called a quotsplinequot The spline is weighted at various points along its length to keep it in contact with the drafting surface as the curve is drawn The simplest implementation is to use linear splines in the intervals between consecutive data points Doing this we denote the approximating function F1x and the quotnquot linear splines which comprise it as f1x This leads to flylx x0 S x S x1 me mx xH s 3c s x 51 nx x 1 S S x Using Newton first order divideddifference interpolating polynomials for the linear splines a1 b1x x0 x0 S x S x1 Fm a woe m xH s c s x 52 an bn x xH x 1 S 5c S xquot The coefficients 61 and b define the individual splines f1 x and are given by a fxH il23n 53 bl fog mi x1 x171 i1 2 3n 53a In the following example linear splines are applied to the data points in Figure 51 Example 51 Suppose y xi in Figure 51 represents the population in millions of a country after xi decades following the initial recorded population Table 51 contains nine population data points yi xi 10000 24570 47709 69495 81951 87151 89031 89675 89988 5 oo cNmJkUJNt O lONUIlkUJNt O O Table 51 Data Points of x in Figure 51 The composite function consisting of linear splines is easily graphed in MATLAB because the plot39 function by default connects data points by linear segments Hence the following MATLAB commands will produce the graph shown in Figure 52 gtgtX07 10 gtgtv0 10 O 10 gtgtFX1 24570 47709 69495 81951 87151 89031 89675 89988 gtgtplcgtt xFx gtgtX1abel x ylabel y gtgtaxisv gtgthcgt1d on gtgtaxis manual gtgtplotxFx The graph in Figure 52 can provide rough estimates of the interpolated values using F1x especially when the zoom39 feature of MATLAB is used However a direct approach utilizes the 39interpl function call which returns a single value or array of linearly interpolated values For example if we wish to estimate the population at the middle of each decade the additional MATLAB commands will do the job man 05 95 gtgtF1nt ntexp1xF7xx1nt xlnt 05000 15000 25000 25000 05000 55000 55000 75000 55000 95000 Frnt17225 3014a 520m 7572 24551 22mm 29253 29727 29231 29930 XEv ya MOPLAX m x t 11 8 1 Karyn Frgure 5 2 Apprornrnatang Funetron for Datarn Table 5 1 Consrstrng ofLmear Sphnes There are srtuataons when the apprornrnatang funetron rs reqmred to pass through eaeh hnear sphnes wru satrsfy the rst requrrernent However unless the data pornts are very taghuy spaeed See Frgure 13 the eornbrned funetaon wru have a zxgrzag look as appearance Ft r L or as they are o en referredto Lrnear splme funcuons are tred down at these knots smoothness refers to the eontanurty ofa funetaon andrts denvataves A funetron rs mfxmtely srnooth at a pornt rf the denvataves of every order are eontrnuous at that pornt dead the rst denvatave ofF1x tn Fr re 5 2 rs ducanunuous at eae knot Unless there are three aohaeent data pornts whteh happen to be eolhnear preeewrse hnear funetrons hke Frrwr11 a1ways enhrbrt a change tn slope at eaeh knot The use of preeewrse second order polyhormals offers a pamal soluaor to the problem xf we restrret the quadratre po1 orma1s to the mt pomts That rs wrth n1 data pomts xx x 012 be 1 quadratre sphrres one per rrrterval you re on enng why eaeh quadratre splme passes through my two eorrseeuave data porrrts even t e could t a quadraae fuhetror aeross subrrrterva1s eompnsed ofthree data porrrts look at Frgure 5 3 2H s 3quot Quadratre 15 2 d Quadratre 1D 5 4quot Quadratae 391 Quadratae 1 2 3 A 5 E 7 Frgure 5 3 Preeewrse Quadraae Polyhormah Through Submtervals ofThree Data Pomts Ft Hr 5 thrrd quadraaes Tab 5 2 hsts the four quadratre po1ynomra1s and therr rst derwaav The MATLAB poly t function was used to obtarr the quadraae polyhormal eoemererrts Interval Quadratre FrrstDenvaave 3x 15x2 05x e5r 235 72 5x2 235x 7 30 73x2 25x 7 28 76x 25 6 2x2 35x 153 4r 5 Table 5 2 Quadraue Fuheuohs m Flgure 5 3 and ThelrFlrst Derlvatlves x2andx4 Frrst Dehyauye of 2rd Quadratle Frrst Derwauye of 1 Ouadraue lsl Derlvahve Frrst Derwatwe of 4quot Quadratle Frrst Derwatwe of 3d Quadraue Flgure 5 4 Frrst Derwatwe of Composrte Fuheuoh Shown rh Flgure 5 3 ose we eohstrarh eaeh quadratre polynomlal to a slngle rhteryal betweeh reasohably s f d t ts the slopes of eaeh rhterror quadratre should be eohtrhuous at both end porhts as well resuluhg rh two more eohstrarhts It would appear that the problem of determrruhg eaeh quadraue spllne 15 now over eohstrarhed by yrrtue freedom represented by the three polynomlal eoemerehts four rhteryals de ned by flve data porhts There are a total of 12 eoemerehts 3 for eaeh of4 qua aue polynomlals to be determrhed hehee there must be xaetly 12 rhdepehdeht eohstrarhts to generate 12 llnearly rhdepehdeht equauohs There are 8 eohstrarruhg of eaeh rnterya1 Eaeh of the 3 rntenor pornts produees on1y a 51 g1e srnoothness eonstratnt equahty of the rst denyatryes of adja nt quadratre sphnes rnereasrng th tota1 nurnber of eonstratnts to 11 Consequenuy the sys tern of hnear equauons for deterrnrnrng the 12 eoemerents 15 not over eonstratned on the eontrary an addmonal ratn re e uauon s needed to rodu rn of equauons wrth a umque so1uuon for the 12 eoemerents The rn 5mg equ uon W111 be reyea1ed m the genera1 dseussron of quadraue sphnes whreh follows a1though you may have gured rt out on your own 1f you looked carefully at the rst quadratre sphne 25 Coef m ents 2U 15 WEI Equahty Constratnts 5 denyatrye El 4 E E 1D 1 Frgure 5 5 A Cornposrte Funetron Fx ofoecewxse Quadraue 5phnes n a set of data pornts xx x 012 h there are n quadraue sphnes to the Newton ddvxdedrddfference form of Equauon 3 7 Smce the srnoothness eonstratnt mvolves the rst denyauye rt 15 not adyrsab1e to ernp1oy the Lagrange form whreh 15 the where the eonstratnts are based so1e1y on the funeuon ya1ues as m Frgure 5 3 where eaeh preeewrse quadratre was requrredto t a subset ofthree data pornts Choosmg the Newton ddvxdedrddfference form the funeuon F20 eornpnsed of quadraue sphnes 15 expressed as shown m Equauon 5 4 Refemng to Frgure 5 o ther are 3 eoemerents to be deterrnrned for eaeh of n sphnes henee a tota1 of 311 eoemerents v s moi aauaxagxqpapqu uomaN mpamasaxdag sauqu axmpenb 9 g am g z 3 x 3H2 KXAXHXAX J quotX4xquotq quotn x 4 042 4 x H2 4 x quotq Hz W S X S X LOCLOWJ Low 4 x s x s 2 240042402 Hm w 0 1 40042 4x z 42 lt1 ix 3 x 3 x ZX4X X4XZJ x4xlq in 2 3 x 3 quotz X4xquot141 g mtg 1 sauqu am no smmsuoa 30 Jaqwnu am e Euxsodmx Aq paumqo ale smaumaoa 358m auxuuaiap 01 511012an The rst set of equations constrain the quadratic splines to the data point function values at the end points of their respective intervals This assures Fzx will be a piecewise continuous function through the entire set of data points I Interval End Point Function Values Must Equal Data Point Values Fzxfx 10 12 n 55 For the rst interval F2x0 a1b1x0 x001x0 x0x0 x1 56 fx0 a 57 Fzx1 a1b1x1 x001x1 x0x1 x1 58 fx1 fxob1x1xo 59 bl M 510 x1 x0 Evaluating the second quadratic spline at the end points of its interval gives Fzx1 a2b2x1 x102x1 x1x1 x2 511 fx a2 512 Fzx2 a2b2x2 x102x2 x1x2 x2 513 fx2 fx1b2x2 x1 514 hi M 515 x2 xt Similar constraints equations are given below F2094 FAX apply to the remaining splines The general form of the a b1 x171 x17101x171 x171x171 x 516 fxH i1 2 3n 517 atQxt x17101xt x171x1 95 518