### Create a StudySoup account

#### Be part of our community, it's free to join!

Already have a StudySoup account? Login here

# INTRODUCTION TO NUMERICAL ANALYSIS MTH 351

OSU

GPA 3.6

### View Full Document

## 17

## 0

## Popular in Course

## Popular in Mathematics (M)

This 125 page Class Notes was uploaded by Miss Johan Jacobson on Monday October 19, 2015. The Class Notes belongs to MTH 351 at Oregon State University taught by Staff in Fall. Since its upload, it has received 17 views. For similar materials see /class/224453/mth-351-oregon-state-university in Mathematics (M) at Oregon State University.

## Reviews for INTRODUCTION TO NUMERICAL ANALYSIS

### What is Karma?

#### Karma is the currency of StudySoup.

#### You can buy or earn more Karma at anytime and redeem it for class notes, study guides, flashcards, and more!

Date Created: 10/19/15

MATLAB NOTES Matlab designed for numerical computing Strongly oriented towards use of arrays one and two di mensional Excellent graphics that are easy to use Powerful interactive facilities and programs can also be written in it It is a procedural language not an object oriented lan guage It has facilities for working with both Fortran and C lan guage programs USING MATLAB At the prompt in Unix type Matab Or use the pull down menu accessed with the right button Run the Demo program simply type demo To seek help on any command simply type help command To seek information on Matlab commands that involve a given word in their description type loo kfor word Look at the various online manuals available thru the course web page MATLAB is an interactive computer language For ex ample to evaluate 3 26 4 72 35 y 96 96 2 use y 6 4X 7XX 3X 5 3X2 There are many built in functions eg eXpX cosX sqrtX logX The default arithmetic used in MATLAB is double preci sion and real However complex arithmetic appears au tomatically when needed sqrt 4 results in an answer of 21 The default output to the screen is to have 4 digits to the right of the decimal point To control the formatting of output to the screen use the command format The default formatting is obtained using format short To obtain the full accuracy available in a number you can use format long The commands format short e format long e will use scientific notation for the output Other format options are also available MATLAB works very efficiently with arrays and many tasks are best done with arrays For example plot sin 90 and cosx on the interval 0 S 90 S 10 t O110 X cost y sint plottxty The statement t azhzb with h gt 0 creates a row vector of the form t aaha2h giving all values a jh that are S b When h is omitted it is assumed to be 1 Thus 11 15 creates the row vector n12345 A R RAYS b 1 2 3 creates a row vector of length 3 A 123456789 creates the square matrix A Mgr 0001M cow Spaces or commas can be used as delimiters in giving the components of an array and a semicolon will separate the various rows of a matriX For a column vector b 13 6J results in the column vector ARRAY OPERATIONS Addition Do componentwise addition A 1 2 3 2 6 1 B 2 3 3 2 2 2 CAB results in the answer 3 5 C 0 0 4 1 Multiplication by a constant Multiply the constant times each component of the array D 2A results in the answer 2 4 D 6 4 Matrix multiplication This has the standard meaning E 1 2 2 1 3 2 F 2 1 3 1 23 G EF results in the answer 1 2 4 5 3 G 2 1i 5 4 3 3 2 8 7 3 A nonstandard notation H 3 F results in the computation 111 2 13 52 H 3i1 1 1i 1 2 3i 2 5 00 I l COM PON ENTWISE OPERATIONS Matlab has com ponent wise operations for multiplication division and exponentiation These three operations are denoted by using a period to precede the usual symbol for the operation With a123 b214 we have ab 2 2 12 ab 05 20 075 a 3 1 8 27 2 a 2 4 8 b a 2 1 64 The eXpression 3 96 2 can be evaluated at all ofthe elements of an array X using y6 4x7x2 35 the command y 6 4X 7XX 3X 5 3X2 The output y is then an array of the same size as X OTHER COMMANDS clear To remove the current variables from use cls To dearthe outputscreen help command name Briefdescription of command name help sqrt resuhsin the output SQRT Square root SQRTX is the square root of the elements of X Complex results are produced if X is not positive Special arrays A zeros23 produces an array with 2 rows and 3 columns with all components set to zero 0 0 0 0 0 0 B ones23 produces an array with 2 rows and 3 columns with all components set to 1 1 1 1 1 1 1 eye3 results in the 3 X 3 identity matrix 1 0 0 0 1 0 0 0 1 ARRAY FUNCTIONS There are many MATLAB commands that operate on arrays we include only a very few here For a vector X row or column of length n we have the following functions maXX 2 maximum component of X minX minimum component of X abs X 2 vector of absolute values of components of X sumX sum of the components of X normX 9612 n2 SCRIPT FILES A list of interactive commands can be stored as a script file For example store t O110 X cost y sint plottXty with the file name pot trigm Then to run the program give the command plottrig The variables used in the script file will be stored locally and parameters given locally are available for use by the script le FUNCTIONS To create a function we proceed similarly but now there are input and output parameters Consider a function for evaluating the polynomial 19013 a1 6021 l a32 anmn l MATLAB does not allow zero subscripts for arrays The following function would be stored under the name polyeval m J ray named coeff and the polynomial is to be evaluated The coefficients a are given to the function in the ar at all of the components of the array X function value polyevalXcoeff function value polyevalXcoeff Evaluate a polynomial at the points given in X The coefficients are to be given in coeff The constant term in the polynomial is coeff1 n lengthcoeff value coeffnonessizeX for i n 1 11 value coeffi Xvalue end l J V V V V V V V V V V V V V NewtonCotes Quadrature Mth 351 Summer 1999 Bent E Petersen Filename newtoncotesmws We investigate various compound closed NewtonCotes quadrature rules I NewtonCotes quadrature proceeds as follows we wish to approximate an integral gt Intfx xa b be x dx 2 We subdivide the interval ab into m the degree subintervals of equal length The corresponding ml endpoints are called the nodes Now let p be the unique interpolation polynomial of degree at most m through the points on the graph of f corresponding to the nodes Then the integral of p is our approximation to the integral of f The Lagrange formula for the interpolation polynomial a change of variable and a quick calculation show that gt Intpx xa b b a Sumfak b a m k0 m Int pk x x0 1 m kb af Z a pkxdx k0 m 0 k l where pk is the interpolation polynomial which is l at and 0 at the other multiples of The m m fpbc dxb a a integrals of the pk are called the NewtonCotes coefficients of degree n Here s a routine to compute these coefficients NCC proc m local khxxyyycoefqw xxseqkmk0m coef yseq0k0m for k from 0 to m do yysubsopk11y make kth component in y 1 in yy qinterpxxyyx wintqx01 coefopcoefw push onto list of weights od coef end Let s check a few Page 1 i Tapezoidal rule 39gt NCC1 Simpson39s rule 39gt NCC2 I Simpson s 38 rule gt NCC3 I Boole s rule 39gt NCC4 F l q 90 45 15 45 90 gt NCC5 19 25 25 25 25 19 288 96 144 144 96 288 39gt NCC6 iii iii 840 35 280 105 280 35 840 39gt NCC7 751 3577 49 2989 2989 49 3577 751 17280 17280 640 17280 17280 640 17280 17280 39gt NCC8 989 2944 464 5248 454 5248 464 2944 989 28350 14175 14175 14175 2835 14175 14175 14175 28350 Something important happened here For the case of degree 8 we have some negative weights When the weights are all positive we can integrate nonnegative functions without worrying about loss of signi cance errors When some weights are negative this is no longer the case As a result the degree 8 method is not popular 39gt NCC9 2857 15741 27 1209 2889 2889 1209 27 15741 2857 89600 89600 2240 5600 44800 44800 5600 2240 89600 89600 gt NCC10 Page 2 16067 26575 16175 5675 4825 17807 4825 5675 16175 26575 16067 598752 149688 199584 12474 11088 24948 11088 12474 199584 149688 598752 Degree 9 is wellbehaved but degree 10 again has a sign problem This problem with the sign of the weights makes it useless to increase the degree in an effort to increase accuracy We may see a decrease in truncation error but an increase in roundoff error Thus compound methods are usually used The idea in the compound methods is we break the interval into n subintervals where n is diVisible by the degree m that we plan to use Then we use the NewtonCotes method on the rst set of m subintervals then on the second set and so on and then add the reuslts together Let39s write a routine to do it for us gt NCCintprocfrrangemn gt local kjabqChhhS gt if typen numeric then gt if not typenposint then ERRORquotn must be positivequot fi gt if not typemposint then ERRORquotm must be positivequot fi gt if not iremnm0 then ERRORquotm degree must divide n orderquot fi gt fi gt aop1r bop2r qnm hban hh baq gt CNCCm S0 gt for j from 0 to q l do gt for k from 0 to m do gt ssCk1fajhhkh gt od od gt Shhs rescale recall intervals here have length b aq gt end 1 For example here is the compound Simpson39s rule with 6 subintervals on ab 39 gt NCCintfa b26 1 1 1 2 5 1 1 2 1 2 1 1 1 1 2 b 61 Ka t a b t a b t a b t a b 3 3 6 3 6 6 3 3 3 3 2 2 3 3 3 2 1 5 1 f ab fb 3 6 6 6 If we multiply by 6 and simplify we see the familiar l424224l pattern Now let39s compare some order 36 methods for gt rng2 5 exIntexpx xrng Page 3 gt V V V 5 ex J ex dx 2 1 36 relerror trapezoidal evalf ex NCCint exp rng 1 36 ex 40 136 relierror trapezoidal 000578636735181ll3819505569654341257633373l 236 relerror Simpson 5 evalf ex NCCint exp rng 2 36 ex 40 236 relierror Simpson39s 2676970493645l23463788605006041036146537 10396 336 relerror Simpson 38 evalf ex NCCint exp rng 3 36 ex 40 336 relierror Simpson 38 60182121005l0396036642660688234988688149 10396 4 36 relerror Boole evalf ex NCCint exp rng 4 36 ex 18 43 6 relierror Boole 70620238582496227l 10399 6 36 relerror evalf ex NCCint exp rng 6 36 ex 40 636 relierror 2470l00308872566719802249620241372804574 103911 9 36 relerror evalf ex NCCint exp rng 9 36 ex 40 936 relierror 2083586305707669970039883286859413078756 103913 12 36 relerror evalf ex NCCint exp rng 12 36 ex 40 1236 relierror 2149309937627081215600849772199922166976 103918 18 36 relerror evalf ex NCCint exp rng 18 36 ex 40 1836 relierror 2967743044315975885121890984893149445570 103925 Note each of these methods required 37 function evaluations In computing the error we use the actual value of the integral and both eX and NCCint are computed symbolically The difference is found and divided by ex also symbolically We then convert to oating point at a prescribed precision By using a very high precision we are able to compute the error with minimal roundoff so we are seeing mostly truncation error here Let s consider another example Page4 gt rng0 Pi exIntsinx xrng ex J sinx dx gt 1 36 relerror trapezoidal evalf ex NCCint sin rng 1 36 ex 40 136 relierror trapezoidal 0006347001875770279795828014456345840185000 236 relerror Simpson 5 evalf ex NCCint sin rng 2 36 ex 40 236 relierror Simpson39s 3224859886427327311172174433400355000000 10396 336 relerror Simpson 38 evalf ex NCCint sin rng 3 36 ex 40 336 relierror Simpson 38 7262524560995178254906258891912465000000 10396 4 36 relerror Boole evalf ex NCCint sin rng 4 36 ex 18 43 6 relierror Boole 938471270000000000 10399 gt 6 36 relerror evalf ex NCCint sin rng 6 36 ex 40 636 relierror 3638559258825456117945759847500000000000 103911 gt 9 36 relerror evalf ex NCCint sin rng 9 36 ex 40 936 relierror 3448427660223488205259756450000000000000 103913 gt 12 36 relerror evalf ex NCCint sin rng 12 36 ex 40 1236 relierror 4440389945798304591145000000000000000000 103918 gt 18 36 relerror evalf ex NCCint sin rng 18 36 ex 40 1836 relierror 8987633535428150000000000000000000000000 103925 V V V One more set of examples gt fx gtexpx rng3 8 exIntfx xrng gt 1 36 relerror trapezoidal evalf ex NCCint f rng 1 36 ex 40 136 relierror trapezoidal 0016069937074547616735101977082522870883l9 gt 236 relerror Simpson 5 evalf ex NCCint f rng 2 36 ex 40 236 relierror Simpson39s 2062533686062336605453489969246517737922 10395 gt 336 relerror Simpson 38 evalf ex NCCint f rng 3 36 ex 40 336 relierror Simpson 38 4630090957512773713123467870871570605425 10395 gt 4 36 relerror Boole evalf ex NCCintfrng4 36 ex18 43 6 relierror Boole 150390032384600922 10397 gt 6 36 relerror evalf ex NCCint f rng 6 36 ex 40 636 relierror l448094774390365865743323428367052296037 10399 gt 9 36 relerror evalf ex NCCint f rng 9 36 ex 40 936 relierror 332592399654838873000877469769565050858l 103911 gt 12 36 relerror evalf ex NCCint f rng 12 36 ex 40 1236 relierror 2568836808982809598089232253667627777154 103915 gt 18 36 relerror evalf ex NCCint f rng 18 36 ex 40 1836 relierror 7011852462380051678651793139847097342379 103921 gt fevalnf unassign For our nice and smooth examples the higher degree methods are much superior The composite Boole39s method is no harder to apply than Simpson39s rule even for hand calculations but has smaller error for smooth lnctions Compound Boole 12 subintervals gt NCCint f a b 4 12 algsubs b a12h 35 190 simplify 90 2 E7fa32fah12fa2h32fa3 h l4fa4h32fa5 h 12fa6h32fa7h14fa8h32fa9h12fa10h32fa11h Page 6 L 7fa12hh Compound Boole 8 subintervals 39 gt NCCintfa b48 algsubsb a8h 190 simplify90 7fa32fah12fa2h32fa3h14fa4h32fa5h 12fa6h32fa7h7fa8hh 2 h We see we have acoef cient of E andthen the sequence 7 32 12 32 14 32 12 32 14 14 32 12 32 7 The 14 s come from two 739s at the points where two of the interpolation intervals meet Since Boole s rule depends on interpolation of degree 4 it is obviously exact for polynomials of degree at most 4 However just as for Simpson39s rule it is actually exact for one degree more that is for polynomials of degree at most 5 gt collectNCCintx gtxquot5a b44 b 1m1J 6 6 The degree 9 rule is one of the higher degree rules which has positive weights We can expect it to be fairly wellbehaved Why don t you experiment a bit with compound degree 9 rule Page 7 Method of Least Squares Mth 351 Oct 10 2001 Maple 6 Bent E Petersen We can solve linear least squares problems by using either of Maple s linear algebra packages but the most convenient approach is to use the stats fit package gt restart gt withstats fitdescribe withstats statplots withplots Warning the name changecoords has been redefined E Suppose we have some experimental data gt Xdata l23456789101112 Xdata z 1 2 3 4 5 6 7 8 9 10 11 12 gt Ydata322314567563 Ydata z 3 2 2 3 1 4 56 7 5 6 8 and we suspect a functional relationship of the form y fX where f contains some adjustable parameters The Maple leastsquare function accepts directly the equation that we are trying to fit to our data Here s a little procedure to compute the L2 norm of a vector or 1 array of arbitrary dimension N2 proc a local k sqrt suma k 2kl nops a end Alternater we could use the norm function from the linear algebra package gt N2 vgtlinalg norm v frobenius E It will also be handy to have a plot of the given data points gt pltO scatterplot XdataYdata colorred thickness3 Constant function fit Consider a constant function gt flxgtc f x gtc 39 gt eqnlleastsquare xy yf1x c XdataYdata l3 eqn 2232 Thus or our optimal functional relationship of the given type is given by gt Flunapplyrhs eqnl x E 3 F12 The statsdescribe package defines the function mean 39 gt ym mean Ydata E ym 3 We see as expected the best least squares fit is the mean of the Ydata The deviation or total residual is the square root of the sum of the squares of the residuals gt NDl N2 Ydatamap Fl Xdata evalf 1 ND 2 474 3 725 71 8035 3 Egt pltlplotFll 12colorblack 7 gt displayplt0pltltitle quotEstimation by the meanquot Estimation by the mean 8 5 4 2 7 o 2 21 6 8 1b 12 The mean is a pretty rough estimator However it is useful to consider since the quantity NDl is used to scale the total residual in other cases to obtain a dimensionless measure of goodness of fit Linear Function Line Fit gt f2xgtmxb f2 22x gt mx b gt ean 1eastsquare xy yf2 x mb XdataYdata 2 74 32 eqquot y 143x 33 gt F2 unapplyrhs ean x 74 32 F2 22x gt x 143 33 gt ND2 N2 Ydatamap F2Xdata evalf 1 ND2 22 2645214 429 3791168734 The dimensionless unadjusted measure of goodness of fit is given by gt gf2lND2N39D1A2 evalf 8214 2 11297 7270956891 Close to l is good close to 0 is bad Our result here is noncomrnital gt p1t2 plot F2 l 12 colorblack 7 gt displayplt0plt2 title quotLeast squares line fitquot Least squares line fit 8 6 4 2 o 9 A 1b 15 Special Line Fit The leastsquare function is quite versatile For example if we have reason to expect the slope should be 12 in the example above then we can restrict our least squares candidates to lines with slope 12 gt f3xgtl2xb 1 f3 x gtExb gt eqn3 1eastsquare xy yf3 x b XdataYdata 3 l 13 e n 2 x 4 y 2 12 gt F3unapplyrhs eqn3 x 1 13 F3 x gt x 2 l2 1 11 11 11 11 11 gt N39D3 N2 Ydatamap F3Xdata evalf 1 ND3 519 6 i 3796928584 gt gf3lN39D3N39D1A2 evalf 3 459 gf 632 7262658228 As expected the fit is not as good as before since we have restricted the class of candidates gt plt3plotF3l 12colorblack 7 gt display plt3plt0 titlequotSpecial line fitquot Special line fit 8 6 4 2 o 2 21 E s 1b 1 7 Quadratic Fit gt f4xgtax 2bxc f4 x gtax2bxc gt eqn4leastsquare xy yf4 x abc XdataYdata 31 2 115 21 X X 1001 1001 11 eqn4 2 y 7 gt F4unapplyrhs eqn4 x i 31 2 115 21 F4 x gt x x 1001 1001 11 39 gt ND4N2 YdatamapF4Xdata evalf 1 ND4 13119106 1001 3618412234 gt gf4l ND4N131 A2 evalf 59420 quot3979079 7514004982 gt plt4plot F4l 12colorblack 39 gt display plt4plt0 titlequotQuadratic fitquot Quadratic fit 8 5 4 2 o 2 21 8 10 1 2 Cubic Fit gt f5xgtax 3bx 2cxd f5x gtax3bx2cxd gt eqn5leastsquare xy yf5 x abcd XdataYdata 56 3 257 2 38251 35 x x x 3861 819 27027 9 eqn5 y I 1 1 VV F5 unapplyrhs ean x quot393003 7 3264388558 gf5l N39D5ND102 evalf V 5 189236 gf quot39237237 7976664686 gt plt5plotF5l 12colorblack V N39DS N2 Ydatamap F5Xdata evalf 1 ND5 1 96098002 display plt5plt0 titlequotCubic fitquot Degree 4 fit gt 6 x gtax 4bx 3 cx 2 dxe Cubic fit 8 5 4 2 o 7 21 E 8 1b 17 f6 x gtax4bx3cx2dxe gt eqn6leastsquare xy yf6 x abcde XdataYdata x x 13728 61776 F6 unapply rhs eqn6 x gt gt N39D6 N2 Ydatamap F6Xdata evalf 61 4 8033 661116 y E 53837 2 X X 41184 61776 198 282169 1319 1 ND6 6476470 858 7 2966073381 7 gt gf6lN39D6N39D1A2 evalf 6 112919 gf quot135564 7 8329571273 gt p1t6 plot F6 1 12colorblack gt display plt6plt0 titlequotDegree 4 fitquot Degree 4 fit 8 6 4 2 o 2 4 8 8 1o 12 Degree 5 Fit gt f7xgtax 5bx 4cx 3dx 2exf f7x gtax5bx4cx3dx2exf gt eqn7 1eastsquare xy yf7 x abcde f XdataYdata 23 5 321 4 82091 3 29543 2 47885 3 x x x x x 10608 4862 116688 9724 9724 11 eqn7 2 y F7 unapplyrhs eqn7 x N39D7 N2 Ydatamap F7Xdata evalf E 1 ND7 y 32344455 gt gt 2339456423 39 gt gf7lN39D7N39D1A2 evalf 344183 384098 8960812084 gt plt7plot F7l 12colorblack gt display plt7plt0 titlequotDegree 5 fitquot Degree5fit r o 2 A s 8 1b 12 Degree 10 Fit V f12 x gta10xquot10a9x 9a8xquot8a7 x 7 a6 xquot6a5x 5a4x 4a3 x 3 a2x 2alxa0 f12 x gta0x10119x9a8x8a7x7a6x6115x5a4x4a3x302x2alx110 eqn12 1eastsquare xy yf12 x a10a9a8a7a6a5a4a3a2al a0 XdataYdata 157 10 14263 9 133939 8 2915323 7 293270623 V eqn12y x x x x x6 3628800 5080320 1693440 2298240 22982400 6533315417 5 84072994307 4 50526430079 3 124616339339 2 278166527 78140160 234420480 x 51279480 x 75969600 x 188955 1055 2 F12 unapply rhs eqn12 x ND12N2 YdatamapF12Xdata evalf 109 y 176358 7786641323 gf12lN39D12N39Dl 2 eva1f 9181259 2 9288188 9884876361 gt plt12plot F12l 12colorblack I l l VV ND12 2 V 7 gt display plt12plt0 titlequotDegree 10 fitquot Degree 10fit 8 6 4 2 7 o 7 1 7 gt This example is typical Even though our goodness of fit is about 0988 the graph intuitively looks worse The fit is of course great near the data points but polynomials of high degree tend to oscillate a lot and therefore we get large unjustified swings between the data points The philosophy more or less is a large goodness of fit is meaningless unless it is achieved by a polynomial of low degree or more generally by a function with few adjustable parameters Is it meaningful in the later case Well maybe maybe not An Exponential Fit Maple can handle any kind of linear least squares fit but not nonlinear gt f20 xgtabexp x15 lle 0xaabe gt eqn20leastsquare xy yf20 x ab XdataYdata gt F20unapplyrhs eqn20 x gt ND20N2 YdatamapF20Xdata evalf 3670066155 gt gf20lN39D20ND1A2 evalf E 7442521727 gt p1t20 plot F20 l 12colorblack gt display p1t20p1t0 titlequotAn exponential fitquot An exponential fit Would you care to use any of these models to predict the quotfuturequot larger values of X What is required to make prediciton anything but a meaningless exercise gt Interpolation Polynomials Mth 351 Spring 2001 Bent E Petersen Filename 351 2001interppoysmws In this worksheet I present a few procedures that are useful in experimenting with interpolation polynomials gt restart f1nterp interpolates a function given the distinct nodes by computing the function values at the nodes and submitting the resulting lists of abscissas and ordinates to Maple interp function Note that Maple allows sending any number of unspeci ed arguments to a procedure so it is a simple matter to pass lists od data of arbitrary length to a procedure All the variables passed to the procedure may be accessed through the list args the nth entry in the list is argsn the sublist consisting of the kth to the nth entry is argskn Note args is not an array in spite of appearances The number of entries in the list of parameters is given by nargs gt finterpprocfx gt if nargs lt 3 then ERRORFAIL fi gt interp args3 nargs mapf args3 nargs x gt end Here s a quick test of f1nterp gt psfinterpsinx0Pi4Pi3Pi22Pi33Pi4Pi gt plot sin x ps x0 Pi axesNORMAL thickness3 titlequot Error in psquot numpoints100 resolution300 Page 1 Evnzrm ps x 1 41615 2 2334725 25 3 4am 7205 4205 M 4905 6905 V fmterpO m m Wm L L ThusI dataxs ccccwcd ms is mef ment of course but allows us to mantam accuracy UnodesO morder Fr um n w m submta39val s gt Unodes rocnab local L k gt y gt if n lt 1 then ERRORFBIL fi gt L gt for k from o to n do LopLakbcan 0d gt opL gt end Wc retum opL m m w a w n v ahstofthenodes Pagzl of taste it allows nargs in other routines to contains the number of nodes plus other parameters and saves us the trouble of counting the number of entries in a formal list Let39s test Unodes gt Unodes6ab 5 l 2 l l l l 2 l 5 a a b a b a b a b a bb 6 6 3 3 2 2 3 3 6 6 Cnodes returns the n1 ChebysheV nodes for an interval ab gt Cnodesprocnab gt local L k gt if n lt 1 then ERRORFAIL fi gt L gt for k from 0 to n do LopLab a1 cos2k1Pi2n22 od gt PL gt end Here s a couple of tests of Cnodes gt Cnodes4 11 vw VS 0 vs W 545 7 gt Cnodes4ab a2ba1 i v5W Jaaba1i4345 ab ailtb agt 1iW ailtb agt1iW J Let s interpolate the sine function on 0 11 with 7 nodes 0dr 7 gt pufinterpsinxUnodes60Pi 8424 x6 25272 x5 x4 x3 3566 x2 308 x 43963 pu 5706 2988 972 5 n6 5 n5 n4 n3 5 n2 5 n 6 5 x5 2916 3294 TC 5 x Bf 1647 3 x2 135 Ex 1728 n4 n3 4 n2 4 11 Page 3 TV thrnr 1 evalfpu16 7 001296680939813 x5 77 01222092994370 r5 7 00641254180168 7 16073505149962 9 7 00284509716638 212 1 00053876843603x pu 7 39 oxEDthickness73titlywarror in gt plotsinx u n points7 7 00 Ermrmw V 0112140508 1 12141618 2 12242528 3 x r 6f 1 M x l demmal p61yr16rrrra1 gt pc s39nxevalfCnodeS60Pi16 x577 0121145395327 006033382849 7 1613673246 x3 int p 77 001285392572 7 002364997394 x2 1000423400 r7 00001331620 We can putthem m standard order by usmg the function 50110 gt sortpCx PagaA 7 001285392572 x5 01211453953 9 7 006033382849 2 1613673246 9 02364997394 x2 1000423400 x7 00001331620 gt plotsinxpc r 39 aoxEDthickness3titlywarror in pcquotnulnpoints100resolutio 00 Ermrin W n 11112000511511 21415152222425281 x smaller than m the umform nodes case gt PagaS COMPUTING ANOMALIES These examples are meant to help motivate the study of machine arithmetic 1 Calculator example Use an HP 15C calculator which contains 10 digits in its display Let 1 2 3 98765 There are keys on the calculator for the mean 1 I L Z 90339 njzl and the standard deviation 8 where 1 n 2 2 8 71 131 3 In our case what should these equal In fact the calculator gives g 98765 8 i 158 VVhy 2 A Fortran program example Consider two pro grams run on a now extinct computer Program A A 10 20 23 B A 10 PRINT AB END Output 10 119E 7 2 23 Program B A 10 20 gtlt 23 SILLY 00 B A 10 PRINT AB END Output 10 00 Why the change since presumably the variable SILLY does not have any connection to B DECIMAL FLOATING POINT NUMBERS Floating point notation is akin to what is called sci entific notation in high school algebra For a nonzero number 90 we can write it in the form 96 Jug106 with e an integer 1 S g lt 10 and a 1 or 1 Thus 50 g 166666 10 101 with a 1 On a decimal computer or calculator we store 90 by instead storing a g and e We must restrict the number of digits in g and the size of the exponent e For example on an HP 15C calculator the number of digits kept in g is 10 and the exponent is restricted to 99 S e S 99 BINARY FLOATING POINT NUMBERS We now do something similar with the binary repre sentation of a number 90 Write as a g 26 with 1 S lt 102 2 and e an integer For example 110 110011001100 2 2 4 a 1 The number 90 is stored in the computer by storing the a g and 6 On all computers there are restrictions on the number of digits in g and the size of e FLOATING POINT NUMBERS When a number 90 outside a computer or calculator is converted into a machine number we denote it by On an HP calculator fl3333 333333333310 10 1 The decimal fraction of infinite length will not fit in the registers of the calculator but the latter 10digit number will fit Some calculators actually carry more digits internally than they allow to be displayed On a binary computer we use a similar notation I will concentrate on a particular form of computer float ing point number that called the IEEE floating point standard In single precision we write such a number as fl a 1a1a2 a232 2e fl a 1a1a2 a232 26 The significand g 1a1a2 a232 immediately sat isfies 1 S g lt 2 What are the limits on e To understand the limits on e and the number of bi nary digits chosen for g we must look roughly at how the number 90 will be stored in the computer Basi cally we store a as a single bit the significand g as 24 bits only 23 need be stored and the exponent fills out 32 bits 4 bytes Thus the exponent e occupies 8 bits including both negative and positive integers Roughly speaking we have that 6 must satisfy 11111112 g e g 11111112 127 g e g 127 In actuality the limits are 126 g e g 127 for reasons related to the storage of O and other num bers such as loo What is the connection ofthe 24 bits in the significand g to the number of decimal digits in the storage of a number 90 into floating point form One way of answering this is to find the integer M for which 1 O lt 90 S M and 90 an integer implies flac 90 and 2 flM17 M1 This integer M is at least as big as 1114 223222320 231 s 2 This sums to 224 1 In addition 224 10 02 224 also stores exactly What about 224 1 It does not store exactly as 24 10 01gt 2 23 0 8 2 Storing this would require 25 bits one more than al lowed Thus M 224 16777216 This means that all 7 digit decimal integers store ex actly along with a few 8 digit integers THE MACHINE EPSILON Let y be the smallest number representable in the ma chine arithmetic that is greater than 1 in the machine The machine epsilon is 77 y 1 It is a widely used measure of the accuracy possible in representing num bers in the machine The number 1 has the simple floating point represen tation 1 1000220 What is the smallest number that is greater than 1 It is 12 23 1001220 gt 1 and the machine epsilon in IEEE single precision float ing point format is 77 2 23 i 119 x 10 7 THE UNIT ROUND Consider the smallest number 6 gt O that is repre sentable in the machine and for which 16gt1 in the arithmetic of the machine For any number 0 lt 04 lt 6 the result of 1 04 is exactly 1 in the machines arithmetic Thus 04 drops off the end39 of the floating point representation in the machine The size of 6 is another way of describing the accuracy attainable in the floating point represen tation of the machine The machine epsilonhas been replacing it in recent years It is not too difficult to derive 6 The number 1 has the simple floating point representation 1 1000220 What is the smallest number which can be added to this without disappearing Certainly we can write 12 231001220gt1 Past this point we need to know whether we are us ing chopped arithmetic or rounded arithmetic We will shortly look at both of these With chopped 2 23 arithmetic 6 and with rounded arithmetic 6 2 24 ROUNDING AND CHOPPING Let us first consider these concepts with decimal arith metic We write a computer floating point number 239 as zaQ106Eaala2an1010 3 with al 75 0 so that there are n decimal digits in the significand a1a2 an10 Given a general number acaa1a2an10106 a1 0 we must shorten it to fit within the computer This is done by either chopping or rounding The floating point chopped version of 90 is given by fl a a1a2 an10 10 3 where we assume that 6 fits within the bounds re quired by the computer or calculator For the rounded version we must decide whether to round up or round down A simplified formula is fl aa1a2an1010 an1 lt 5 a a1a2 an10 00110 10 an1 2 5 The term 00110 denotes 10n1 giving the or dinary sense of rounding with which you are familiar In the single case 00 0an1an2 10 00 0500 10 a more elaborate procedure is used so as to assure an unbiased rounding CHOPPINGROUNDING IN BINARY Let 201a2an22 with all a7 equal to O or 1 Then for a Chopped floating point representation we have fl a 16L2 an2 26 For a rounded floating point representation we have fl 01a2an2106 an1 0 0 1a2 an2 OO12 10 3 an1 1 ERRORS The error x flx 0 when 90 needs no change to be put into the computer or calculator Of more interest is the case when the error is nonzero Consider first the case 90 gt 0 meaning a 1 The case with 90 lt O is the same except for the sign being opposite With 90 75 and using chopping we have fl lt 96 and the error 90 flac is always positive This later has major consequences in extended numerical com putations With 90 75 flac and rounding the error 90 flac is negative for half the values of ac and it is positive for the other half of possible values of ac We often write the relative error as x fl x This can be expanded to obtain W lt1 an Thus flx can be considered as a perturbed value of 90 This is used in many analyses of the effects of chopping and rounding errors in numerical computa tions For bounds on a we have 2 5 2 rounding 2n1 5 0 chopping IEEE ARITHMETIC We are only giving the minimal characteristics of IEEE arithmetic There are many options available on the types of arithmetic and the choppingrounding The default arithmetic uses rounding Single precision arithmetic 7124 126 e 127 This results in M 224 16777216 77 2 23 119 x 10 7 Double precision arithmetic n 53 1022 g e g 1023 What are M and 77 There is also an extended representation having 71 69 digits in its significand NUMERICAL PRECISION IN MATLAB MATLAB can be used to generate the binary float ing point representation of a number Execute the command format hex This will cause all subsequent numerical output to the screen to be given in hexadecimal format base 16 For example listing the number 7 results in an output of 4010000000000000 The 16 hexadecimal digits are 0 1 2 3 4 5 6 7 8 9 abcdef To obtain the binary representation convert each hexadecimal digit to a four digit binary number For the above number we obtain the binary expansion 0100 0000 0001 1100 0000 0000 for the number 7 in IEEE double precision floating point format NUMERICAL PRECISION IN FORTRAN In Fortran variables take on default types if no explicit typing is given If a variable begins with I J K L M or N then the default type is INTEGER Otherwise the default type is REAL or SINGLE PRECISION We have other variable types including DOUBLE PRE CISION Redefining the default typing Use the statement IMPLICIT DOUBLE PRECISIONA HO Z to change the original default of REAL to DOUBLE PRECISION You can always override the default typ ing with explicit typing For example DOUBLE PRECISION INTEGRAL MEAN INTEGER P Q TMLOUT FORTRAN CONSTANTS If you want to have a constant be DOUBLE PRECI SION you should make a habit of ending it with DO For example consider DOUBLE PRECISION Pl Pl314159265358979 This will be compiled in way you did not intend The number will be rounded to single precision length and then stored in a constant table At run time it will be retrieved zeros will be appended to extend it to double precision length and then it will be stored in Pl lnstead write Pl314159265358979D0 Simpson39s Rule and Cubics Mth 351 July 14 1999 Bent E Petersen Filename simpsoncubicmws Degree 3 and 2 It is well known that Simpson39s quadrature is more precise than one at rst expects and moreover is exact for cubics The second fact accounts for the rst one We can use Maple to demonstrate symbolically that Simpson39s rule is exact for cubics Morever we can demonstrate a similar phenomenon for higher degree polynomials and other compound quadrature rules based on interpolation To show that Simpson39s rule is exact for cubics it suffices to show that for any three points on a cubic with equispaced abscissas the area under the graph of the cubic from the rst to the last point is the same as the area under the graph on the same interval of the unique quadratic throught the three points We begin by constructing a general cubic polynomial lnction gt p3 unapply axquot3bxquot2cxd x p3 x gt ax3bxzcxd Now let39s compute the integral from r h to r h for any real numbers r and h gt A3 collect int p3 x xr h rh h 2 A3 b2arh32ar32d2br220rh Now we compute the quadratic interpolation polynomial p2 through the points on the graph of p3 with abscissas r h r and r h gt X2r hrrh X2 r hrrh ar h3br hzcr hdar3brzcrdarh3brh2crhd gt q2unapplyinterp X2Y2x x q2x gt3arbxzah2 3ar2cxar3 arh2d Let39s compute the integral from r hto r hfor any real numbers r and h gt B2collectintq2 x xr h rh h 2 BZ b2arh323arbr22 3ar2cr2ar32dh gt simplify A3 B2 0 It is now easy to see that Simpson39s compound quadrature rule that is Simpson s rule is exact for cubics i i i gt Y2mapP3X2 i i i i l Page 1 Degree 5 and 4 Let s verify that similar results hold for higher degree We rst show that given any quintic polynomial degree 5 and ve points on its graph with equispaced abscissas the area under the graph of the quintic from the rst to the last point is the same as the area under the graph over the same interval of the unique quartic degree 4 interpolation polynomial through the given points We begin by constructing a general quintic polynomial gt p5 unapply axquot5bxquot4cxquot3dxquot2exfx p5 x gt ax5bx4cx3dx2exf Now let39s compute the integral from r 2 h to r 2 h for any real numbers r and h gt A5 collect int p5 x xr 2h r2h h 64 5 160 3 l6 2 3 145 b64ar h Tard32br l6cr h 4ar54cr34br44drz4er4fh Now we compute the quartic interpolation polynomial p4 through the points on the graph of p5 with abscissas r 2hr h r rhandr2h gt X4r 2hr hrrhr2h X4r 2hr hrrhr2h Il Il ll ll Il II II I gt Y4mapp5X4 Y4ar 2h5br 2h4cr 2h3dr 2h2er 2h ar h5br h4cr h3dr hzer h ar5br4cr3drzer arh5brh4crh3drhzerh ar2h5br2h4cr2h3dr2hzer2hf gt q4unapplyinterp X4Y4x x q4x gtb5arx4c 10ar25ah2x3 15arh210ar3dx2 4ah4elSarzhz 5ar4x 5ar3hzfar54arh4 Let s compute the integral from r 2 htor2 hfor any real numbers rand h gt B4 collect int q4 x xr 2h r2h h 64 5 2 160 3 l6 2 3 34 b64ar h l6c lOarrTard32b5arr h 4b5arr44ar5410ar3dr24c 10ar2r34e 5ar4r4fh gt simplify A5 B4 0 The algebra would be tedious and error prone to do by hand but Maple has no problem with it so let s try one more example Degree 7 and 6 Page 2 Now let39s verify that similar results hold for degree 7 and 6 We show that given any polynomial of degree 7 and seven points on its graph with equispaced abscissas the area under the graph of the degree 7 polynomial from the rst to the last point is the same as the area under the graph over the same interval of the unique interpolation polynomial of degree 6 through the given points We begin by constructing the general polynomial of degree 7 gt gt p7 unapply axquot7bxquot6cxquot5dxquot4exquot3fxquot2gxquot1h x p7x gt ax7bx6cx5dx4ex3fx2gxh Now let39s compute the integral from r 3 h to r 3 h for any real numbers r and h gt A7 collect int p7 x xr 3h r3h h 74 43 7 3 2 486 5 A7 4374ar b h 3402ar l458br d486cr h 54er 180cr3378ar5 108 er 18f270 br4h36hz 6ar76er36fr26br66gr6dr46cr5h Now we compute the interpolation polynomial p6 of degree 6 through the points on the graph of p7with abscissas r 3 h r 2hr h r rh r2handr3 h gt X6r 3hr 2hr hrrhr2hr3h X6r 3hr Zhr hrrhr2hr3h gt Y6mapp7X6 Y6ar 3h7br 3h6cr 3h5dr 3h4er 3h3fr 3h2 gr 3hhar 2h7br 2h6cr 2h5dr 2h4er 2h3 fr 2h2gr 2hh ar h7br h6cr h5dr h4er h3fr hzgr hh ar7br6cr5dr4er3frzgrh arh7brh6crh5drh4erh3frhzgrhh ar2h7br2h6cr2h5dr2h4er2h3fr2h2gr2h har3h7br3h6cr3h5dr3h4er3h3fr3h2 gr3hh gt q6unapplyinterp X6Y6x x Il Il Il ll Il II Il ll q6x gt7arbx6 21ar2cl4ahzx5al 70arhz35ar3x4 35ar4e 49ah4l40arzh2x3f21ar5 140ar3h2l47arh4x2 36ah6g70ar4h2 l47arzh4 7ar6xh49ar3h4 36arh6ar7 l4ar5h2 Let39s compute the integral from r 3 hto r 3 h gt B6collectintq6 x xr 3h r3h h 4374 7 36 4374ar 7 b h Page3 486 486 21ar2cr3402ar3dl4587arbr2h5378ar5 180 21ar2cr32707arbr4108d35 ar3r254 35 ar4er18 h3 6h267arbr66ar76f21ar5rz6 21ar2cr56 7ar6gr 6d35ar3r46 35ar4er3h gt simplify A7 B6 0 Are you beginning to suspect a general fact Note the direct approach taken above is not suitable for hand calculation even for fairly low degree Eventually Maple will give up too But no nite number of cases will suf ce to prove the general fact anyway For a proof what we really need is a suitable error estimate for the compound interpolation quadrature rules with equispaced nodes Page 4 EVALUATING A POLYNOMIAL Consider having a polynomial p a0 a1a22 ann which you need to evaluate for many values of 90 How do you evaluate it This may seem a strange question but the answer is not as obvious as you might think The standard way written in a loose algorithmic for mat poly a0 for 1 n poly poly ajacj end To compare the costs of different numerical meth ods we do an operations count and then we compare these for the competing methods Above the counts are as follows additions n 1 multiplications 1 2 3 n This assumes each term ajxj is computed indepen dently of the remaining terms in the polynomial Next do the terms acj recursively 963296907 1 Then to compute 2903 will cost n 1 mul tiplications Our algorithm becomes poly a0 alas power 2 96 for 2 n power 2 ac power poly poly aj p0wer end The total operations cost is additions n multiplications n n 1 2n 1 When n is evenly moderately large this is much less than for the first method of evaluating For ex ample with n 20 the first method has 210 multi plications whereas the second has 39 multiplications We now considered nested multiplication As exam ples of particular degrees write n 2 p a0 a1 agx n 3 PW a0W1 a2a3 n 4 PW a0W1 a2a3a4 These contain respectively 2 3 and 4 multiplica tions This is less than the preceding method which would have need 3 5 and 7 multiplications respec tively For the general case write PW a0W1 a2Wn 1 anquot39 This requires n multiplications which is only about half that for the preceding method For an algorithm write poly 2 an forjn 1 1O poly aj 96 p0ly end With all three methods the number of additions is n but the number of multiplications can be dramatically different for large values of n NESTED MULTIPLICATION Imagine we are evaluating the polynomial p a0 a1a22 ann at a point 90 239 Thus with nested multiplication 19Z a0za1Za2Zan 1anz 39 We can write this as the following sequence of oper ations bn 2 an bn l an l l39 25b an 2 l39 an l w 3 b0 2 a0 l zbl The quantities bn1 b0 are simply the quantities in parentheses starting from the inner most and working outward Introduce q 2 b1 13296 133962 moan 1 Claim 1996 190 96 ZCI Proof Simply expand 130 96 239 131 3296 133962 moan 1 and use the fact that ijij1 aj1 j1n With this result we have P 96 50 we 90 z 90 239 Thus qx is the quotient when dividing 1990 by x z and b0 is the remainder lf 239 is a zero of then b0 2 O and then 1996 96 zq For the remaining roots of we can concentrate on finding those of In rootfinding for polynomi als this process of reducing the size of the problem is called deflation Another consequence of is the following Form the derivative of with respect to ac obtaining Mm 96 zq 6196 p z C1Z Thus to evaluate 1990 and px simultaneously at 90 239 we can use nested multiplication for 192 and we can use the intermediate steps of this to also evaluate pz This is useful when doing rootfinding problems for polynomials by means of Newton s method APPROXIMATING SF Define 1 CI sin 75 SFac dt 90 75 0 96 0 t We use Taylor polynomials to approximate this func tion to obtain a way to compute it with accuracy and simplicity As an example begin with the degree 3 Taylor ap proximation to sin 75 expanded about t O 1 1 sint t t3 t5 cosct 6 120 with ct between 0 and t Then sint 1 1 1 t2 t4 t 6 120 cosct 13 39 t 513 1 1 2d7 1 t2 t4cosct dt 0 t 0 6 120 1 3 1 t4 dt COSC 18 120 0 t 133sintdt 1 1 2R 96 m 0 t 18 2 1 1 COS Ct 0 How large is the error in the approximation 1 SF 1 1 82 on the interval 11 Since cosct S 1 we have for 90 gt O that 1 1 t4dt 12096 0 1 4 600 and the same result can be shown for 90 lt 0 Then 0 S R2 for S 1 we have 0 lt R lt 1 x 2 600 To obtain a more accurate approximation we can pro ceed exactly as above but simply use a higher degree approximation to sin 75 In the book we consider finding a Taylor polynomial approximation to SF with its error satisfying lR8l s 5 gtlt 10 9 Incl 1 A Matlab program plotsintm implementing this approximation is given in the text and in the class account The one in the class account includes the needed additional functions sinttaym and polyevenm Begin with a Taylor series for sin 75 3 5 39tt 1n 539quot 35 H 2n 1 t27 Li 1 1n 2n1cosct with ct between 0 and t Then write 1 m t2 t4 S39t 1 quot3996 m 0i 35 1 t2n 2 1 71 dt R WAN 2n 2m x2 x4 2 Q m 1 2n 2 1 R a 2n R2n2 1nm cosct dt 1 CB t2n R22x g 0 1 m cosctdt To simplify matters let 90 gt 0 Since cosct S 1 96271 R lt 1 m 15 d t l 2 2l 0 2n l 1 2n12n1 It is easy to see that this bound is also valid for 90 lt 0 As required choose the degree so that IR2n2xl 5 gtlt 10 9 From the error bound 1 max R 90 lt lmlgl l 2 2 l 2n 12n 1 Choose n so that this upper bound is itself bounded by 5 X 10 9 This is true if 2n 1 Z 11 Le n 2 5 The polynomial is 2 4 6 8 1 1lt lt1 pm 3355 7799 JC and ISFx pxl s 5 gtlt 10 9 Incl 1 To evaluate it efficiently we set u 2 and evaluate 1 u u2 u3 u4 9W 1 8 E 35280 3265920 After the evaluation of the coefficients done once the total number of arithmetic evaluations is 4 addi tions and 5 multiplications to evaluate 1990 for each value of ac Mth 351 Least Squares Mth 351 Summer 2002 July 28 2002 Maple 5 and 6 Bent E Petersen Filename 351u2002leastsquaresmws gt restart Maple has support for least squares fitting in the stats fit package and also in the linalg package The stats fit package is the more convenient one to use for fitting data to a given equation gt withstats fit It will also be useful to have some nice plot commands available Warning the name changecoords has been redefined We will need some random data to illustrate the commands E gt withplots gt Xdata seqk2k0 20 1 3 5 7 9 11 13 15 17 19 Xdata 01112131415 6 7 8 9 10 2 2 2 2 2 2 2 2 2 2 gt Ydata map xgtxsin x evalf Xdata Ydata z 0 9794255386 1841470985 2497494987 2909297427 3098472144 3141120008 3149216772 3243197505 3522469882 4041075725 4794459674 5720584502 6715119988 7656986599 8437999977 8989358247 9298487113 9412118485 9424848880 9455978889 These data points all he on the graph of X sinX Here s a plot of the graph for comparison with the least squares fits that we obtain below gt plot xsin x x0 10 tit1equotxsinxquot thickness3colorblue xsinx Let s try a leastsquares line fit gt eqnl ymxb eqnlymxb gt lsql leastsquare xy eqnl mb XdataYdata lsql 2 y 9789783571x 2636407548 Note Xy tells Maple the names of the variables in eqnl and mb tells Maple which parameters to adjust If you want to leave b as an undetermined parameter you could do it as foolws gt lsqlbleastsquare xy eqnl m XdataYdata 13 qu 2 y l4634l4634 b 101755993lx b To plot the solution lsql we note we just want to plot the right hand side of the equation lsql Let s plot it together with X sinX so we can compare them gt plot xsin x rhs lsql x0 10color bluered thickness3 tit lequotlsq line and xsinx II Isq line and xsinx Of course it makes more sense to compare our least squares fit with the data points First we need the Cartesian coordinates corresponding to our data gt L U for k from 1 to nops Xdata do L op L Xdata k Ydata k M ed 7 gt evalfL4 0 0 5000 9794 1 1841 1500 2497 2 2909 2500 3098 3 3141 3500 3149 4 3243 4500 3522 5 4041 5500 4794 6 5721 6500 6715 7 7657 7500 8438 8 8989 8500 9298 9 9412 9500 9425 10 9456 imgO PLOT POINTS op L SYMBOL CROSS 20 COLOR RGB 0 0 1 imgl plot rhs lsql x0 10 colorred thickness3 titlequotLine fitquot display img0 imgl VV ll Il II V Line fit Looking at the data you may feel a polynomial of degree 4 ought to fit better Let s try it gt eqn2 yabxcxquot2 dx 3 ex 4 651112yabxcxzdx3ex4 gt lsq2zleastsquare xy eqn2 abcde XdataYdata lsq2 2 y 07174541678 2973013639 x 1138367815 x2 1940921577 x3 01008454155 x4 gt img2 plot rhs lsq2 x0 10 colorred thickness3 titlequotPolynomi al of degree 4quot 7 gt display img0img2 Polynomial of degree 4 Don t confuse correlation and causation The degree 4 fit certainly looks better but you should be careful not to draw any unwarranted conclusions Extrapolation is definitely out For example the degree 4 polynomial here decreases after x 10 whereas x sinx continues to increase If you were dealing with experimental data the good fit might mislead you You really need to have an underlying model which implies a certain functional relation before a least squares fit with that function allows you to reasonably extrapolate or draw any conclusions other than statistical ones concerning properties of the actual data gt Elementary Numerical Analysis Contents 1 TAYLOR Polynomials to Errors Lo Numeric Differentiation 1 Roots Midpoint Method u Roots NEWTON s Method ax Roots HALLEY s Method q Roots Secant Method 00 Roots Miscellaneous w Fixed Point Iteration 0 NEWTON sMethod for Systems 1 Secant Method for Systems 2 Interpolation polynomials Bent E Petersen Summer 2004 7 20040806 13 Interpolation polynomials and NEWTON divided differences 19 a Interpolation and EBYsEV polynomials 21 15 Interpolation Splines 22 16 Least Squares Fitting 24 17 Compound NewtonCotes Quadrature M ethods 00 RICHARDSON Extrapolation and ROMBERG Qudrature 32 w GAUSSLEGENDRE Quadrature 34 to 3 Miscellaneous Quadrature Problems 34 to Iterative methods for xed points for non linear systems Simultaneous and succes sive displacements to to Vectors Inner Products and Norms 38 to Lo M atrices 39 B E Petersen Version 20040806 Elementary Numerical Analysis 24 GAUSSJORDAN Row Reduction 41 26 Positive De nite Matrices 43 27 GAUSSJACOBI and GAUSSSEIDEL iter 25 LU Facrorization 42 ative methods for linear systems This list of problems will be augmented during the quarter Some of the problems are from old lists of sample problems old tests or old assignments Thus there may be some duplication or at least near duplication of other material on my web pages Other problems may seem familiar because I used them as examples in my lectures Some of the problems here may be used with minor editing on future tests Other problems espe cially those which are pretty theoretical are for instruction and amusement If you have dif culties with any of the problems then ask me for hints or clari cation but don t wait until the end of the term 1 TAYLOR Polynomials Problem 11 We know from calculus that EULER s number e is given by e M8 51H k H o Deduce25 S e S 275 Hint Show k 2 23k 2 ifk 2 3 Problem 12 Suppose we use the TAYLOR sum 1 27 Me 51quot a H 0 to approximate E elQ How large should we take n to guarantee the error is no larger than 10 Hint You may nd it useful to note e12 3 2 Problem 13 Suppose we approximate 52 on the interval 71 1 by the the TAYLOR polyno mia 4 px1 295 2x2 E953 Use the TAYLOR remainder to give a good upper bound for the truncation error that we make at x 75 an at x 7 Problem 14 Suppose we approximate x logz 1 on the interval 71 1 by the the TAYLOR polynomial px x 7 952 953 Mm 351 2 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Use the TAYLOR remainder to give a good upper bound for the absolute value of the truncation error l x 7 that we make at x 7 and at x Note 1097 is a correct upper bound but it is not a good one Problem 15 Let px be the TAYLORiMACLAURIN polynomial of degree 5 for the function x sinx expx A Compute B Suppose we use px to estimate x on the interval 0 1 Use TAYLOR s remainder to bound the error in p1 and in Use the estimate 5 S 3 in computing your error bounds Problem 16 We approximate a certain function x by an nth degree TAYLOR polynomial about the origin If the absolute value of the kth derivative of x on the interval 0 1 is bounded by 4 C for each k how large whould we take n in order to estimate 05 with an error no larger that 02 in absolute value Problem 17 Suppose O 0 f 0 2 f 0 7 and lf xl S 0024for0 S x S 2 Estimate 1 by using a TAYLOR polynomial of degree 2 Compute a good bound for the absolute error Problem 18 Suppose O 1 flt1gt0 72 flt2gt 0 1 flt3gt 0 71 and lf4xl g 0384 for S 1 A Estimate 05 by using a TAYLOR polynomial of degree 3 3 B Compute a good bound for the absolute error Problem 19 Suppose l 1 flt1gt1 72f2gt1 1 flt3gt1 71 and flt4gt1 24 Estimate 15 by using a TAYLOR polynomial of degree 3 4 Problem 110 Find the TAYLOR polynomial with center at 0 and degree 4 for the function f given by flx ec s 1 xz Given that l f5gt l 3 0005 for l xi 3 05 estimate the errorin approximating 05 by means of its TAYLOR polynomial of degree 4 Compare your estimate with the actual error Problem 111 Suppose 0 1 f 0 72 f 0 7 f 0 g and lf4xl 3 0016 for 0 S x S 2 Estimate 1 by using a TAYLOR polynomial of degree 3 Compute a good bound for the absolute error Problem 112 Suppose 1 1 f 1 1 f 1 72 f 1 4 Estimate 125 by using a TAYLOR polynomial of degree 3 Problem 113 We approximate a certain function x by an nth degree TAYLOR polynomial about the origin If the absolute valuse of the kth derivative of x on the interval 71 1 is bounded by 2 C for each k how large should we take n in order to estimate 025 with an error no larger than 000002 Problem 114 We know 3 3 5 x x x xiggsinx xii forxgt0 i Mth 351 3 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Use this inequality and only this inequality to obtain an estimate of 1 sinx dx 0 x With an error no bigger than 00008 Be sure to explain What you are doing and to back up your claim that the error is no bigger than 00008 Hint Find an interval containing the value of the integral and use the midpoint as your estimate Problem 115 Part A Let x em The TAYLOR polynomial of degree 6 for x is 1 1 1 1 1 727747i5776 x 2x 8x 159 2409 Given that fm S 120 for S estimate the error in approximating x by the TAYLOR polynomial on the interval 7 Part B Compute the value of the TAYLOR polynomial at x Part C Use your calculator to compute f and compare this value With the estimate obtained in part B How does the actual error compare With your estimate of the error obtained in part A Problem 116 Compute the TAYLOR polynomial of degree 4 With center at 1 for Evaluate the TAYLOR polynomial at x 12 and compare your result With V12 as given by a calculator Compare the actual error With the error estimate given by the TAYLOR remainder Problem 117 Substitutet 732 in 1 71tt2t3 t4 1717 lit and integrate With respect to s to obtain x3 x5 967 as 33 arctanxx7i7O H7013 Now use tan 7 1 to estimate 7r Note that the absolute error in your estimate of 7r is boundedby 11 3 8 S 2 d3 0 1 3 How does this estimate compare With the actual error 3 1 8 8 33 d3 3 gwii 19 3 0003190556 0 Problem 118 Estimate 7r as in the previous problem but this time use 1 t7 71 t t2 t3 t4 t5 t6 lit 17t Again compare the estimated error and the actual error Mm 351 4 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Problem 119 Show n 2 2quot11fn 2 2 Then use to deduce e S 3 Problem 120 Recall quot 1 1 e5 Z 7961c 7 egxn1 l l k0 k n 1 for some f between 0 and x If we use the polynomial part here the Taylor polynomial of degree n to estimate e how large must we take n in order to ensure the absolute truncation error is bounded by 10 9 on the interval 71 1 How large must we take n in order to ensure the same error bound 1 1 but just on the interval 7 5 5 Problem 121 If we have a full TAYLOR series we can use it to estimate the error in a Tay lor polynomial This method is sometimes more convenient the estimating the TAYLOR remainder directly and avoids the mysterius f For example if we integrate 220 ark we obtain 7l0g17x Zixk lt 1 lC1 Thus n 00 1 1 flog17 x Z Eyck R7IZI where R7IZI Z 7x k1 kn1 An easy estimate yields 1 Rnw i x k 1 lle n1k 1 n117lxl Taking x we obtain quot 1 log2 2 k2 1c1 a 1 with an absolute error bounded by n How large should we take n to guarantee an estimate of log2 with an absolute error no larger than 1 12 Problem 122 Replacing 95 by 795 in the previous problem yields 00 71k1 log1x2xk lxllt1 lC1 Thus k 1x 00 952 1 loglt17xgt 7log17xlog1x22ki lxllt1 k0 Mm 351 5 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Takingxweobtain 2 1 1 1 27 7 7 0g 3Z2k1lt9gt Estimate the absolute error l Rn l in a H o where has an absolute error no larger than 10 Problem 123 The Taylor polynomial of degree 6 with center at the origin for the function exsinx is given by 7 l 7 37 12 227374775779 pm a x 6x 39 120 36096 One can show 1 7 fquot x g 0095 for 0 g x g 06 71 Use this fact to estimate the error in approximating f05 by p05 Simplify The actual error is actually quite a bit smaller Ans 00007421875000 Problem 124 For a certain function f we have the TAYLOR polynomial of degree 7 at the origin is given by 7 2 2 1 3 10 4 px7 x i x t x Suppose we know 72 s We 3 0 for 71 S x S 1 Use the TAYLOR POLYNOMIAL to estimate f12 Then nd a bound for the error that we make in approximating f12 by p12 Problem 125 Let x log1 x and note that k 7 1 k 71H kgt1 f z lt gt HM Thus the TAYLOR polynomial of degree 8 with center at the origin is x2 x3 x4 x5 x6 x7 3 ptzgtxiiiiiiiiiii 2 3 4 5 6 7 8 Find an upper bound for the absolute error in px as an approximation to log1 x for 0 S x S by using the TAYLOR remainder or facts about alternating series C Use the TAYLOR remainder to obtain an upper bound for the absolute error inpx for 7 S x S Mm 351 6 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis 2 Errors A real number you is said to have m signi cant decimal digits as an approximation to the real number act y 0 provided that w lt 5 X10 m 1 9 The best value of m is called the number ofsigni cant decimal digits in you as an approximation to act or simply the number of signi cant digits when act and decimal are understood It is clearly given by m flog10 lt2 Here l j indicates the greatest integer function otherwise known as oor Note this de nition of signi cant decimal digits may not agree with the count obtained by other methods but it is more convenient and equally useful xt xa Problem 21 The fractions 272 and 5T are wellknown approximations to 7r Use the formula above to nd the the number of signi cant decimal digits in each Does your understanding of signi cant digits agree with your calculated result Problem 22 Imagine an impractical binary oating point representation with exponent rang ing from 763 to 64 and with a mantissa of length 12 bits Assume only normalized representations are used and there is no packing A What is the largest exactly representable positive number B What is the smallest exactly representable positive number C What is the unit round Answers A 253 212 7 1 m 369 X 1019 B 2 63 m 108 X 10 19 C 2 11 m 488 X 10 4 Problem 23 In the previous problem suppose we implement gradual under ow by allowing denormals when the exponentis 763 Suppose also we reserve the largest exponent 64 to designate too projective in nity NaN s and other special data How does this affect the answers to the previous problem Problem 24 Imagine an impractical binary oating point representation with exponent rang ing from 763 to 64 and with a mantissa of length 12 bits Assume only normalized representations are used and the format is packed that is the most signi cant bit in the mantissa is not explicitly stored Mth 351 7 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis A What is the largest exactly representable positive number B What is the smallest exactly representable positive number C What is the unit round Problem 25 In the previous problem suppose we implement gradual under ow by allowing denormals when the exponentis 763 Suppose also we reserve the largest exponent 64 to designate too projective in nity NaN s and other special data How does this affect the answers to the previous problem Note that numbers with exponent 763 are not packed to allow for denormals Problem 26 Imagine an impractical binary oating point representation with exponent rang ing from 763 to 64 and with a mantissa of length 12 bits Assume only normalized representations are used and there is no packing Find the error that we make when storing 1132 decimal Ans 180 00125 Problem 27 In the previous problem if the length of the mantissa is only 7 bits what would be the error in storing 1132 3 Numeric Differentiation Problem 31 Here is an unusual second order method for approximating the rst derivative fa w 0012 Show that the error term can be written as cflt3gt an 9 I13 and compute the constant c Hint Use Taylor expansions Problem 32 Here is an unusual rst order method for approximating the second derivative fa fa73fa2h22fa3h 3h Show that the error term can be written as cflt3gt 11 9012 and compute the constatnt c Hint Use TAYLOR expansions Here are a few standard numeric differentiation formula Mm 351 8 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis FIRST DERIVATIVE M h 7 fta H OW W 9012 central 3point fa W H 9012 forward 3p0mt 2h i a 2h 8fa high 8fa 7 h u 7 2h 9014 central 5point S ECOND DERIVATIVE f a W 9012 central 3point f a i a 3h 4fa 7 5fa h 2fa 9012 forward 4point THIRD DERIVATIVE fHa u 2h 7 2101 hygf a 7 h 7 u 7 2h 9012 Central 5Point W 7 73fa4h14fa3h724fa2h18fah75fa f a 7 2m 0012 forward 5point Problem 33 Given rs 234 f19 221 f20 212 andf21 207 estimate f 2 0 and f 2 0 by using central second order estimates Problem 34 The function f has the values shown in the table below H 95 HH 0 H 0125 H 0250 H 0375 H 0500 H 0625 H0750 H 0875 H 1000 H H x HH 2802 H 2873 H 2911 H 3109 H 3002 H 2999 H 2981 H 2786 H 2544 H Estimate f 0 875 and f 0 875 by using second order central estimates Problem 35 The function f has the values shown in the table below H 95 HH 0 H 0125 H 0250 H 0375 H 0500 H 0625 H0750 H 0875 H 1000 H H x HH 2802 H 2873 H 2911 H 3109 H 3002 H 2999 H 2981 H 2786 H 2544 H What estimate of f 0500 do you get from the fourth order central 5point method listed above with step size 0125 Use the same step size and compute the estimate given by order 2 central 3point method Explain why there is no reason to expect one result to be better than the other for this problem Mm 351 9 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis 4 Roots Midpoint Method Problem 41 Show analytically that the polynomial pxx376x29x75 has a root in the interval 4 5 Suppose we bisect the interval A Which half of the interval can we guarantee contains a root B If we use the midpoint of the subinterval containing the root to estimate the root give an upper bound for the error Problem 42 The polynomial px x373x1 has a root in the interval 0 1 Use the bisection method also called binary search method to locate the root with an error no larger than Your work must show clearly that you used the bisection method Problem 43 Letpx 2x3 7 3x2 6x 5 The polynomial px has a rootin the interval 71 0 A Suppose we bisect the interval Which half of the interval can we guarantee contains a root Why B Bisect the subinterval known to contain a root Which subsubinterval can we now guarantee contains the root C If we use the midpoint of this last interval as an estimate of the root nd an bound for the error that we make Problem 44 Show analytically that the polynomial px 2x4 7 495 1 has a root in the interval 0 1 Suppose we bisect the interval A Which half of the interval can we guarantee contains a root B If we use the midpoint of the subinterval containing the root to estimate the root give an upper bound for the error Problem 45 Show analytically that the polynomial px x5 x 1 has a rootin the interval 710 705 Suppose we bisect the interval A Which half of the interval can we guarantee contains a root B If we use the midpoint of the subinterval containing the root to estimate the root give an upper bound for the error Problem 46 Let x x cosx Part A Show that the function f has exactly one root in the interval in2 0 Part B Use the bisection method to obtain an estimate of the root with an error no larger than 7r8 In part B you must give the estimate provided by application of the bisection method and explain it I have no interest in the more accurate solution provided by your calculator 5 Roots NEWTON s Method This section deals with the NEWTON root estimate iteration fxn 9 I may Mth 351 10 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Problem 51 The polynomial px x5 x 1 has a root in the interval 71 0 Use the midpoint of this interval as an initial guess and apply NEWTON s method once What is your new approximation to the root Problem 52 The polynomial px x373x1 has a root in the interval 0 1 Use the midpoint of this interval as an initial guess and use one iteration of NEWTON s method to obtain a new approximation to the root You must obtain exactly the value given by one iteration of NEWTON s method Problem 53 Letpx x3 x2 3x 7 4 Let x0 1 be an initial guess to a root A Use N EWTON s method twice to compute successive approximations x1 x to a root B Estimate the error in the root estimate 951 C Given that px has only one real root 866369759 t t i nd the actual errors in 951 and 952 How do they compare with your estimate of the error in 951 Problem 54 Find two interval of length 1 with integer endpoints each interval guaranteed to contain a root of at expx 7 695 Be sure to justify your answer without nding the roots Use an initial guess we 0 and apply NEWTON s iteration once to obtain a new estimate x1 for one of the roots Problem 55 Let a gt 0 and let x x3 7 1 Suppose we decide to estimate the cube root 113 by applying NEWTON s method to estimate the positive root of For n 2 0 nd an expression for the n 1 iterate yen in terms of 957 Simplify Problem 56 One of the roots of the polynomial px x2 7 x 7 1 is the golden ratio 1 16180339 t t t Use NEWTON s method with initial guess we 2 for the root and compute the iterates x1 x2 and x3 and also the error in each iterate Is the rate of convergence about what you would expect Problem 57 Let px x5 7 3x3 4x 5 Let x0 72 be an initial guess to a root Use NEWTON s method to compute successive approximations x1 x2 x3 x4 and x5 Given that 71 627772774213681160 t t t is the only real root ofpx nd the actual errors in the root estimates Do they behave as you would expect Problem 58 Let x7x6 74x5 76x42x3 8x25x1 and let I g x 7 l f tx so 9 has the same roots as f but the roots of g are all simple Use an initial guess we 1 1 and apply NEWTON s method to f and to 9 Comment on the rate of convergence for example compare the errors to the root i 1 6180339887498948482045868343t t t a double root for f Mth 351 11 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Problem 59 Let x7x6 74x5 76x42x3 8952 5x1 The polynomial f has a root of multiplicity 3 at x 71 Take an initail guess of 950 709 and apply NEWTON s method and the modi ed NEWTON s method xmxquot 3g3 and compare the errors In each case comment if the convergence appears to be linear or quadratic Problem 510 Use NEWTON s method to approximate one of the two roots of 1 2 ftx10g96 gm 10gI x Take for your initial guess we 50 and iterate three times to obtain x3 In applying NEWTON s method we may not wish to compute the derivative at each step Thus we could try freezing it say fxn 96 x T W This idea is probably not a good one in the scalar case here but it has advantages in dealing with solving nonlinear systems Problem 511 The polynomial px x3 7 3x2 2x 7 4 has a simple root at approximately 27963219032594415351 Use NEWTON s method with an initial guess we 30 and 5 iterations to estimate the root Compute the error in each iterate and comment on the convergence Problem 512 Repeat the previous problem but using the modi ed NEWTON s method de scribed above with the derivative frozen at 30 Compare the apparent rates of convergence 6 Roots HALLEY s Method In class we derived NEWTON s method from the tangent line approximation to a graph We can also derive it directly from a TAYLOR expansion If 95 is an approximation to a root and we write the root as 95 6 we have 0 ftxn 5 ftxn f n5 for some 7 between 95 and 95 6 Thus 6 7 m Ma and we expect when all goes well 6 m 67 where fxn 67 7 f Mm 351 12 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Accordingly the next root estimate is ftxn xn1xn6nxn7 HALLEY s method results from taking a higher order TAYLOR expansion 1 0 fm 6 7 me mm Ef TnnW form some 717 between 95 and 95 6 This equation is quadratic in 6 We could use the quadratic formula but we would probably get a headache Instead note the 6 which satis es f 957 6 0 is given by f M Ma where of course we do not know 5 If we substitute this expression for one of the factors of 6 in the 62 term we obtain o m mm 7 j g l It follows that 2ftxnf n meow5 7 jawm Thus when all works well we expect 6 m An where 7 2mm 7 ftxnf xn Setting 95 95 An as the next root estimate we obtain HALLEY s iteration m mums 1 xn1 7 In an 1 21W If we start suf ciently close to a simple root 1 then some heavy algebra and analysis shows HAL LEY s method yields An 3 Ka 7 xn1 m Ka 7 for a certain constant K that is HALLEY s method is of order 3 In more prosaic terms if we start suf ciently close to a simple root then after things settle down we will see roughly a tripling of the number of signi cant digits with each iteration Problem 61 The increasing function f de ned by we 7 x 7 case has a simple root given approximately by a 0739085133215160641655312087673873404013411758900757465 Starting with an initial guess of 950 01 compute a few NEWTON and HALLEY iterations and the error in each Comment on your results Are your results consistent with your expectations Note The error in 954 for NEWTON s method is about 7103 X 10 11 and for HALLEY s method it is about 237 X 1049 If you are using the limited precision of an ordinary calculator there is not much point in going past 952 or 953 in this example Mm 351 13 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis 7 Roots Secant Method Problem 71 The polynomial px x373x1 has a root in the interval 0 1 Use the points 950 0 0 and x1 05 as initial guesses and use one iteration of the secant method to obtain a new approximation to the root You must obtain exactly the value given by one iteration of the secant method Problem 72 Letpx x5 7 3953 495 5 Let x0 7210 and x1 7118 beinitial guesses to a root Use the secant method to compute successive approximations x2 x3 x4 and x5 Given that 7116277727742136811601 t t is the only real root of px nd the actual errors in the root estimates Do they behave as you would expect 8 Roots Miscellaneous Problem 81 The polynomial px x5 7 7x4 10952 7 1000 has a simple root at x 5 Thus for small 6 gt 0 the polynomial px 6954 has a root at x 16 which is continuously differentiable in 6 and satis es 10 5 Compute the derivative 1 0 Problem 82 The polynomial px x5 7 7x4 10952 1000 has a simple root at x 5 Thus for small 6 gt 0 the polynomial px 6954 has a root at x 16 which is continuously differentiable in 6 and satis es 10 5 Compute the derivative 1 0 Problem 83 Let fx expx 7 x2 Prove that f has precisely one real root Explain why the root must be in the interval 71 0 Hint Sign change Problem 84 The polynomial px x3 x2 7 x 7 1 has a simple root at x 1 Thus for small 6 gt 0 the polynomial px 6952 has a root at x 16 which is continuously differentiable in 6 and satis es 10 1 Compute the derivative 1 Use your result to estimate a root of px x3 711952 7 at 71 9 Fixed Point Iteration Problem 91 Let 1 1 7 996 95 Show that g has a unique positive xed point 1 Let x0 2 Does the iteration 95 gxn appear to converge to 1 Problem 92 IfGx 95 then G maps theinterval I 9 49 into itselfand l G x l S g for each x E I Part A How do you know that G has a unique xed pointin 9 49 Part B If we 49 and we de ne 96n1 0067 Mth 351 14 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis can you guarantee that the iterates xn converge to the xed point Part C Algebra Find the exact xed point Problem 93 Plot y cossinz and y z Estimate the point where the graphs intersect as an estimate of a xed point for cossinz Deduce also that cossinz has a unique xed point by showing as 7 cossinz is increasing Do a few xed point iterations to improve your xed point estimate Remark the xed point is about 0768169156 i i Problem 94 Plot y coscosx and y z Estimate the point where the graphs intersect as an estimate of a xed point for coscosx Deduce also that coscosx has a unique xed point by showing 957 coscosx is increasing Do a few xed point iterations to improve your xed point estimate Remark the xed point is about 0739085133 i Problem 95 The function f given by x 95 xl has a xed point at 1 and another xed point near 2 Graphically estimate the xed point near 2 What happens if you run the xed point iteration with initial guesses 20 and 19 Estimate the xed point by applying a root nding iteration to approximate a root of x 7 z Ans 1 93456321 i i 10 NEWTON s Method for Systems Warning The problems in this section are very heavy computationally Do not plan on zipping through them The frozen coef cient matrix method referred to below is just NEWTON s method with the partial derivatives evaluated just at the initial point and not updated in subsequent iterations Problem 101 Use NEWTON s method for systems to approximate a solution to the system of equations 952 7 2x 7 y 77 952 432 4 Use the initial point 20 0 3 and iterate twice Given that the actual solution is 1900676726 03112185654 compute the error in each of your iterates Comment What happens if you use the initial point 7100 025 Do three iterations Problem 102 Repeat the previous problem but using the frozen coef cient matrix algo rithm Do several more iterations Comments Problem 103 Use NEWTON s method for systems to approximate a solution to the system of equations 33 7 2952 395 74 Mm 351 15 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Use the initial point 70804 and iterate three times Given that the actual solution is 70865034729 0450902472 compute the error in each of your iterates Comment Problem 104 Repeat the previous problem but using the frozen coef cient matrix algo rithm Do several more iterations Comments Problem 105 Use NEWTON s method for systems to approximate a solution to the system of equations x2y223 3 x2y722 0 x3y37322 2 Use the initial point 10 7100 and iterate three times Given that the actual solution is 144825 7095663 7023250 compute the error in each of your iterates Comment Problem 106 Repeat the previous problem but using the frozen coef cient matrix algo rithm Do several more iterations Comments Problem 107 The solution of the system x2y2 1 952 71 is obviously 0 71 Use the initial point 025 720 and iterate 5 times At each step compute the error in your iterate What behavior of the errors do you observe Problem 108 Repeat the previous problem but using the frozen coef cient matrix algo rithm Do several more iterations Comments Problem 109 Estimate the minimum point of x4 2y4 3952312 29533 3x 7 2y2 4 by setting the partial derivatives equal to 0 and using NEWTON s method With 3 iterations and an initial guess of 710 1 0 Problem 1010 Repeat the previous problem but using the frozen coef cient matrix algo rithm Do several more iterations Comments Problem 1011 The system of equations 2x2 y2 4 x2 7 495g 7 y2 0 Mm 351 16 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis has four solutions Use NEWTON s method with the initial point 10 1 0 and with four iterations to estimate one solution What happens if you use the initial points 710 1 0 1 0 710 or 710 71 0 Problem 1012 Repeat the previous problem but using the frozen coef cient matrix algo rithm Do several more iterations Comments 11 Secant Method for Systems The secant iteration for approximating roots can be extended to systems of equations This method has the advantage of not requiring derivatives but it does require us to supply 3 initial points in the case of 2 variables To estimate a solution to the system aw 0 9062 0 we rst choose 3 points 951 yl x2 yg and x3 yg which are close to a solution By evaluating f at these 3 points we obtain 3 points in 3space xn yn x yn n 1 2 3 If these points are not collinear there is a unique plane P1 through them Likewise from 9 we obtain 3 points 957 yn 9xn yn n 1 2 3 If these points are not collinear there is a unique plane P2 through them If the planes P1 and P2 are not parallel they intersect in a line and if this line is not parallel to the x yplane it meets that plane in a point x4 y4 Now use the points 952 yg x3 yg and x4 y4 and repeat the construction to obtain the next iterate x5 y5 Problem 111 Use the secant method for systems to approximate a solution to the system of equations 953 7 3y 2 0 y372x23x4 0 Use the initial points 708 04 708 05 and 709 05 and iterate three times Given that the actual solution is 70865034729 0450902472 compute the error in each of your iterates Com ment Remark As a check on your work note the rst iterate is 708646900900 04492075016 12 Interpolation polynomials Problem 121 The function fx satis es f0 1 f1 1and f2 2 A Find the interpolation polynomial of degree at most 2 which interpolates f between these points B Use the result of part A to estimate C What estimate do you get for if you use linear interpolation between successive points Mm 351 1 7 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Problem 122 The function x satis es 71 1 O 2 and 2 1 A Find the interpolation polynomial of degree at most 2 which passes through the given points on the graph of B Use this interpolation polynomial to estimate 1 C If it is known that the third derivative of x on the interval 71 2 is bounded in absolute value by 112 nd a bound for the error in your estimate of 1 Problem 123 Compute the interpolation polynomial of degree at most 2 for the function x x3 with nodes at 71 0 and 1 Compute also the interpolation polynomial in the case the nodes are 0 1 and 2 Problem 124 The function x satis es O 1 1 land 2 2 A Find the interpolation polynomial P2x of degree at most 2 which interpolates f between these points B Suppose we know l gl S 20 for 0 S t S 2 Compute a good bound for the absolute error in P205 as an approximation to 05 Problem 125 The function x satis es O 1 1 0 2 0 and 3 71 A Find the interpolation polynomial P3x of degree at most 3 which interpolates f between these points B Suppose we know 1f 3 1502 for 0 S t S 3 Compute a good bound for the absolute error in P315 as an approximation to 15 Problem 126 Compute the interpolation polynomial of degree at most 2 for the function x x3 with nodes at 71 0 and 1 Compute also the interpolation polynomial in the case the nodes are 0 1 and 2 Problem 127 The function x satis es O 1 1 1and 2 2 A Find the interpolation polynomial px of degree at most 2 which interpolates f between these points B Suppose we know l gl S 20 for 0 S t S 2 Compute a good bound for the absolute error in p05 as an approximation to 05 Problem 128 The function x satis es O 1 1 0 2 0 and 3 71 A Find the interpolation polynomial px of degree at most 3 which interpolates x between these points B Suppose we know lf4tl 3 1502 for 0 S t S 3 Compute a good bound for the absolute error in p15 as an approximation to 15 Problem 129 The function x satis es O 3 1 1 and 2 2 A Find the interpolation polynomial of degree at most 2 which interpolates f between these points B Use the result of part A to estimate C What estimate do you get for if you use linear interpolation between successive points Problem 1210 The function x satis es 71 1 O 2 and 2 1 A Find the interpolation polynomial of degree at most 2 which passes through the given points on the graph of B Use this interpolation polynomial to estimate 1 C If it is known that the third derivative of x on the interval 71 2 is bounded in absolute value by 112 nd a bound for the error in your estimate of 1 Problem 1211 Find the interpolation polynomial px of degree 3 3 through the points Lil 23 32 4r Mth 351 18 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Problem 1212 Find theinterpolation polynomial of degree 3 3 through the points 0 2 1 0 2 3 and 3 1 Ans 72x3 175952 7 x 2 Problem 1213 Given the following data compute the interpolation polynomial of degree 3 5 with the given nodes As a check on your work compare your result with my graph below Interpolation Polynomial 13 Interpolation polynomials and NEWTON divided differences DIVIDED DIFFEREN CES flfo il 961 7x0 awth jack flx1quot39axkxllifi0quot39xkill H06 f960fl960a961l9696039fl960961 xnlhi olm i nill Pn71xfl960x1quot39 xnlx fol 39 39 96 Mir where Pk is the interpolation polynomial of degree 3 k for the function f with nodes at 960961quot 39 axle Mth 351 19 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis A remarkable identity which holds for any function f is 1 7 Pntx flxo 961 m 96196 7 MW 7 961 x 7 In A consequence of this identity is if we view Pnx as an approximation to x then the error is given by 1 Wfn1 95 950x95 951 39 39 x i 957 for somefbetween minx0 x1 xn x and maxx0 x1 xn Problem 131 Let x0 0 x1 1 5 and x2 2 and suppose xo 1 x1 0 and x 71 Compute the NEWTON divided differences fx0 x1 and fx0 x1x2 Find the linear interpolation polynomial P1 through the rst two points and the quadratic inter polation polynomial P2x through all three points What estimates do these polynomials give for f0 Problem 132 Let x x3 Let x0 x1 and x be three distinct points Compute the NEWTON divided difference fx0 x1 x2 Problem 133 For a certain function x we know the NEWTON divided differences f7 1 2 f71 1 1 f71 1 2 72 f71 1 2 4 2 A Find the interpolation polynomial P3x for x of degree at most 3 with nodes at 71 1 2 4 B Use the interpolation polynomial to estimate 0 C Find a good upper bound for the absolute error in the estimate of 0 given that f71 1 2 4 3 14 for all x in the interval 71 4 Problem 134 Let x anx 111x7 1 a1xa0 be a polynomial of degree n and let x0 x1 xn be any n 1 distinct points Use properties of the NEWTON divided differences to compute fx0 x1 xn Problem 135 Carefully ll in the missing entries in the following table of NEWTON divided differences Use the NEWTON divided differences to compute the interpolation polynomial of degree at most 3 for the4 datapoints xk yk k 0 3 Problem 136 Fill in the three missing entries in the following table of NEWTON divided differences Mm 351 20 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Use the NEWTON divided differences to compute the interpolation polynomial of degree at most 3 for the4 datapoints zk yk k 0 3 Problem 137 Let x0 0 x1 15 x 2 and x3 3 and suppose fx0 1 ag 1 an 15 and ags 4 Compute the NEWTON divided differences aco aco x1 aco x1 x2 and aco x1 x2 953 Then nd the cubic interpolation polynomial px for x with nodes x0 x1 x2 and x3 Problem 138 Let x 953 Let x1 x2 x3 and 954 be four distinct points Compute the divided difference fx1x2x3x4 Problem 139 For a certain function f we have the NEWTON divided differences flll 2afl152l 3 fl1a2a4l 2 fl1a2a4a5l 5 flla2a4a395a6l 4 Findtheinterpolationpolynomialpx ofdegreeg 4throughthepoints k fk k 1 2 4 5 6 14 Interpolation and EBY EV polynomials the nth EBYsEV polynomial Tnx is given by T x cosnt9 ifx 6089 One can show T006 T x x T2x 2x2 7 1 T3x 4x3 7 3x T4x 8x4 7 8952 1 One can show Tn1x 2x T x 7 Tn1x Let 214T be the normalized manic EB Y EV polynomial Notice N 7 1771 mngnxl72 The roots of Tnx all lie in the interval 71 1 The nth EBYsEV nodes in 71 1 950 951 xn are the roots of Tn1x They are given by 2k 1 xk7cosltm7rgt k01n Mm 351 21 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis For other intervals we dilate and translate appropriately Thus for the interval 1 b the EBYsEV nodes are given by ab ail 2k1 xkT 2 coslt2n27rgt k01 Problem 141 Let 1 x 160z2 Construct theinterpolation polynomial for x with 13 equispaced nodes in theinterval 71 1 and with 13 SEE Y EV nodes so n 12 relative to the interval 71 1 View each interpolation poly nomial as an approximation of f and estimate the maximum absolute error in 71 1 graphically Ans As a check on your work here is the interpolation polynomial with equispaced nodes 922640625 12 4695215625 10 34832075625 3 7326034875 6 7 I 465796 I 931592 9 7452736 9 3726368 2848365225 4 59056665 21 7452736 9 18631849 If you use a calculator your result may contain some odd degree terms with small coef cients due to roundoff The error in the equispaced case is about 80 In the CEB Y EV case the error is much smaller 15 Interpolation Splines Problem 151 Suppose x3ax274xg 03 39 7 3 2 796 996 bx34 2 Find constants a b c such the 896 is twice continuously differentiable on the interval 0 4 Problem 152 Let x0 0 x1 1 5 962 2 and x3 2 5 and suppose fx0 1 ag 1 an 1 5 and ags 1 Computethe NEWTON divideddifferencesfx0 aco x1 aco x1 x2 and fx0x1x2x3 Then nd the cubic interpolation polynomial P3 for x with nodes x0 x1 962 an 963 Problem 153 Suppose 2 x3bx 03 Find all values of the constants a I such the 896 is twice continuously differentiable on the interval 72 2 Problem 154 If 896 0 for x lt 2 and 896 x 7 23 for x 2 2 is it true that s is a cubic spline Justify your answer Mm 351 22 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Problem 155 Let lixax2x3 if0 x l S 3bxcx27x3 iflltx 2 Determine a b and c so that 595 is a natural cubic spline on the interval 0 2 Problem 156 Let x3ax2 if0 x l 596 2 xibxl 1flltx 2 Determine a and b so that 595 is a cubic spline on the interval 0 2 Problem 157 Let lixax2x3 if0 x l S 3bxcx27x3 iflltx 2 Determine a b and c so that 595 is a natural cubic spline on the interval 0 2 Problem 158 Deleted Problem 159 Let 0 11 fo sac x733 if2ltx ls 595 a cubic spline Justify your answer Problem 1510 Deleted Problem 1511 Let n7 x3ax2 if0 x l 7 whim if1ltz 2 Determine a and b so that 595 is a cubic spline on the interval 0 2 Problem 1512 A Write a careful de nition of a cubic spline With nodes at 951 lt x lt lt z B Let 7 6x7x13 ifo 71 m7 6xx13 if 71ltx ls 595 a cubic spline Justify your answer Problem 1513 Let 595 be the natural cubic interpolation spline With knots 0 2 1 0 2 3 and 3 1 Suppose 271731xgx3 if0 zlt1 395 7xl5x27x3 ifl xlt2 733 34xAx2Bx3 29533 Mth 351 23 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis FindAand B Ans A 715 and B Problem 1514 Given the following data do a natural cubic spline t to the data Note the calculations can be done by hand but are rather lengthy You can do the problem easily in Maple but of course you have to nd out how rst Have fun As a check on your work compare your result with my graph below Interpolation Spline 16 Least Squares Fitting Problem 161 Consider a vehicle travelling at 2 mph that is v fps Applying the brakes develops a coef cient of friction a If the vehicle stops in a distance of 3 feet the work done by friction is nmgs where m is the mass of the vehicle and g is the acceleration of gravity say 32 ftsecQ This work must be equal to the initial kinetic energy m 12 Thus if the speed is in mph then the stopping distance in feet is given by 7 121 v2 5 3600 a If the driver s recognition of the need to stop and reaction time to apply the brakes is t then a distance of vt is travelled before the brakes are even applied It follows that the total stopping distance is 22 121 1 2 77v 5Ev 3600M Mm 351 24 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis If the road surface is good the brakes and tires are in good shape and stopping is achieved without locking the brakes a coef cient of friction M of 07 to 08 may be achieved The recognitionireaction time t is probably around 1 second perhaps 2 or more seconds if the driver is tired Suppose we determine the following stopping distances experimentally for a given driver and vehi cle Use the method of least squares to t this data to the relation 3 av bv2 Once you have computed a and I use these values to determine the recognitionireaction time t and the coef cient of friction M Problem 162 A thermistor is a device which can be used as a temperature sensor since its resistance to an electrical current changes with temperature It is necessary to rst calibrate it For a thermistor the theoretical relation between the absolute temperature in Kelvins T that is the temperature in Celcius plus 27315 and the resistance in Ohms R is given by RAe or 1 T ablogR where log R is the natural logarithm for certain constants A and B or a and b In practice the relation is not so tidy A popular relation used in practice is the empirical SteinhartHart equation a blogR c 00ng where c is usually quite small compared to a and b To calibrate a thermistor means to determine the coef cients a b and c experimentally by measuring the resistance at various temperatures and then tting the SteinhartHart relation While the method of least squares would be reasonable often just 3 data points are used and the corresponding linear equations are solved The calibration results are usually provided by the manufacturer who will probably also warn that the SteinhartHart relation is unreliable for temperatures outside the range used in calibration a sure sign of ad hoc tting For examples of real SteinhartHart coef cients see http WWW atpserisor comritcsteirihartsteinhart html Now both R and T have a distinguished 0 point so we do not need to worry about translation invari ance but it is hard to believe that Ohms and degrees Kelvin have some special physical signi cance In other words the form of the equation ought to be invariant under a change of units A change of units for example Ohms Q to kiloOhms k9 in the SteinhartHart relation will introduce a log R2 term On the principle that the form of the relation should be independent of the choice of Mth 351 25 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis units we could introduce this extra term from the beginning Thus we have the modi ed Steinhart Hart relation 1 T c1 CQ log R 005513 2 C4log R3 In general nothing is really gained by this re nement since the SteinhartHart relation is simply a conventional deeply ingrained way to interpolate approximately between measured data points and gives good results in practice Given the following experimental data for a given thermistor calibrate it that is use least squares tting to determine the SteinhartHart coef cients a b and c m 1 in Kelvins M 32919 31482 30741 30210 t 29815 29502 28837 Resistancein Ohms H 8000 15000 20000 25000 30000 35000 50000 Note the resistance at 25 C 29815quot K in the above case 30 k9 is called the reference resistance of the thermistor In the next ve drill problems you are asked to compare the values at the nodes of a least squares approximation with the given data It would be best to do the comparison graphically if possible Problem 163 Given the following data In I do a least squares t to the equation y a bx Compare the values of this least squares polynomial at the nodes with the given values Problem 164 Given the following data In I do a least squares t to the equation y a bx c952 Compare the values of this least squares polynomial at the nodes with the given values Problem 165 Given the following data Mm 351 26 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis In I do a least squares t to the equation yabxcx2dx3 Compare the values of this least squares polynomial at the nodes with the given values Problem 166 Given the following data In I do a least squares t to the equation y ago bags cx5 Compare the values of this least squares polynomial at the nodes with the given values Problem 167 Given the following data In do a least squares t to the equation yabxcx2 l de ex4fx5 Compare the values of this least squares polynomial at the nodes with the given values What hap pened Why Problem 168 Suppose we wish to t the data to the equation y M eAx by using the method of least squares We nd that we need to solve an unpleasant system of nonlinear equations To obtain a more tractable problem we transform the equation to logy logM A95 and solve the corresponding linear least squares problem don t forget to transform the data too Mm 351 27 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis This method is known as transformation It is important to realize that it is an act of desperation The solution of the transformed problem in general will not agree with the solution of the original problem assuming we can even nd it Solve the transformed problem exponentiate to obtain the desired relation and compare the values at the nodes with the given original data Problem 169 Given the following data In I do a least squares t to the equation 1 a bx by rst transforming the problem to a linear least squares problem Compare the values of the least squares relation at the nodes with the given values graphically if possible y Problem 1610 Given the following data In I do a least squares t to the equation 33 a bx c952 by rst transforming the problem to a linear least squares problem Compare the values provided by the least squares equation at the nodes with the given values graphically if possible As a check on your work compare your result with my graph below Cube root of least squares quadratic Mth 351 28 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Problem 1611 Find a and b so that the relation y ago b logx best ts the data points 1 2 2 4 3 7 in the sense of least squares Ans a 195201 b 078422 17 Compound NewtonCotes Quadrature Methods NEWTONCOTES numeric quadrature to approximate x dot of degree m with n subintervals where n must be divisible by m may be described as follows Let the stepsize h be given by b 7 a h 7 n and let one akh k01n Now let Pj be the interpolation polynomial of degree 3 m with nodes 071mz fxj1mla 239 0 1 39 quot am forj 1 Then the compoundNEWTONCOTES approximation is given by 7bm mm N001 m f Z PM dz j1 071M Note when m n the word compoundis dropped If the error is 9011 then we say that the method has order p If the method is exact for polynomials of degree q then we say the method is degree q exact From the construction it is clear that NCn m f is degree m exact but in fact if m is even then the method is degree m 1 exact It follows then that NCn m f is of order m 1 if m is odd and of order m 2 if m is even NCn 1 f is called the compound trapezoidal rule or the compound 2point rule We have N001 Lf 3 lf960 2fx1 2fx2 2fxn2 2fxn71 fxn NCn 2 f n even is called the compound SIMPSON s rule or the comound 3point rule We have New 2 f 2 mm 4fx1 2M 4M 2M4 41 to Mth 351 29 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis NC n 3 f n divisible by 3 is called the compound SIMPSON s 38 s rule or the compound4point rule We have NC 3 f 3541 lflxol 3f951 3f952 2f953 3f964 3f955 3fxn71 f95nl NCn 4 f n divisible by 4 is called the compound BOOLE s rule or the compound 5point rule We have 2h NCn 4 E 7fx0 32fx112fx2 32fx314fx4 32fxn71 MM NCn 5 f n divisible by 5 is called the compound 6point rule We have N001 5 f 19fx0 75f961 50f962 501 751 381 75f966 75fxn1 Since the formulae above after the rst few become somewhat unwieldy it is preferable to simply give the sequence of coef cients for the noncompound method Thus h trapezoidal rule 1 1 SIMPSON s rule 14 1 3h SIMPSON s 38 s rule 1 3 3 1 2h Booms rule E732 1232 7 a 288 h 7point m41216272722721641 7h 8point WW1 3577 1323 2989 2989 1323 3577 751 4h 9point m989 5888 7928 10496 74540 10496 7928 5888 989 6point 19 7550 50 75 19 Note the minus signs in the 9point formula Problem 171 We know 3 e doc m 19 0855369231876677409 u 0 Mm 351 30 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Estimate the integral using the trapezoidal SIMPSON s SIMPSON s 38 s Boole and 7point rules each with 12 subintervals That is compute each of the estimates NC121f NC122f NC123f NC124 f and NC126f where x e on the interval 0 3 Then com pute the error in each estimate Comments Remark As a check on your work note the error in NC12 6 f is 72888341 X 10 7 and the other errors are larger in absolute value Problem 172 Use SIMPSON s rule with 5 function evaluations 4 intervals to estimate 7r f0 s1nx clx F1nd the actual error Problem 173 Use the trapezoidal rule with 5 function evaluations 4 intervals to estimate 0 sinx clx Find the actual error Problem 174 A Use SIMPSON s rule with 2 subintervals 3 function evaluations to es timate 02 x4 clx B Use SIMPSON s rule with 4 subintervals to estimate the same integral C 1 RICHARDSON s extrapolation suggests that R4 0694 7 S2 may be a better estimate Com pute R4 D Compute the actual error in each of the above parts 1 cl 4 x2 01x Use SIMPSON s rule with 4 subintervals to obtain a rational approximation to 7r ie an approxima tion of the form 7r m2550 Estimate the error Problem 175 We know Problem 176 Use SIMPSON s rule with 5 function evaluations 4 intervals to estimate 0 sinx clx Find the actual error Problem 177 Use the trapezoidal rule with 5 function evaluations 4 intervals to estimate 0 sinx clx Find the actual error Problem 178 The function f has the values shown in the table below t x H 0 t 0125 t 0250 t 0375 t 0500 t 0625 0750 t 0875 t 1000 t i x H 2802 t 2873 t 2911 t 3109 t 3002 t 2999 t 2981 t 2786 t 2544 t Estimate the integral fol x clx as follows A Using the trapezoidal rule with 4 subintervals B Using Simpson s rule with 2 subintervals C Using Simpson s rule with 4 subintervals In all four cases round your answer to 3 decimal places Problem 179 Let Tn be the compound trapezoidal rule with n subintervals so n1 nodes Let Sn be the compound SIMPSON s rule with n subintervals n even Suppose for a certain function f on an interval 1 b we have T3f 8812138 T6f 5893109 and T120 5109739 Compute SS f and Su Ans 4920099 4848616 Problem 1710 Use T2 to estimate ff x4 dx Ans 57 Problem 1711 Use S2 to estimate 24 x5 dx Ans 676 Mm 351 31 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis 18 RICHARDSON Extrapolation and ROMBERG Qudrature Problem 181 Consider the elementary integral 1 x4dx 02 0 Let T1 T2 and T4 be the trapezoidal estimates of this integral for 12 and 4 subintervals respectively If we apply RICHARDSON s extrapolation we obtain the SIMPSON s rule estimates S2 and S4 If we extrapolate once again we obtain BOOLE s rule B4 3 Compute T1 T2 T4 S2 S4 and B4 and the corresponding actual errors b Comment on the observed errors Problem 182 Consider the integral 7r2 60895 clx 20000 42 Let T2 and T4 be the trapezoidal estimates of this integral for 2 and 4 subintervals respectively If we apply RICHARDSON s extrapolation we obtain the SIMPSON s rule estimate S4 Compute T2 T4 and S4 Problem 183 Use SIMPSON s rule with 3 function evaluations to obtain an estimate S2 of the integral 0 sinx clx Use SIMPSON s rule with 5 function evaluations to obtain an estimate S4 of the integral 07f sinx clx RICHARDSON s extrapolation suggests that R4 16S4 7 S2 may be a better estimate than S4 Compute R4 and the error in each of S4 and R4 Problem 184 A Use SIMPSON s rule with 2 subintervals 3 function evaluations to es timate 02 955 also B Use SIMPSON s rule with 4 subintervals to estimate the same integral C RICHARDSON s extrapolation suggests that R4 16S4 7 S2 may be a better estimate Com pute R4 D Compute the actual error in each of the above parts E Do your results agree with your expectation Problem 185 Consider the elementary integral 1 x4dx 02 0 Let T1 T2 and T4 be the ccompound trapezoidal estimates of this integral for 12 and 4 subinter vals respectively If we apply Richardson s extrapolation we obtain the Simpson s rule estimates S2 and S4 If we extrapolate once again we obtain BOOLE s rule B4 A Compute T1 T2 T4 S2 S4 and B4 and the corresponding actual errors B Comment on the observed errors Problem 186 Consider the integral 7r2 60895 clx 20000 42 Let T2 and T4 be the compound trapezoidal estimates of this integral for 2 and 4 subintervals respectively If we apply Richardson s extrapolation we obtain the SIMPSON s rule estimate S4 Compute T2 T4 and S4 Mm 351 32 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Problem 187 Use SIMPSON s rule with 3 function evaluations to obtain an estimate S2 of the integral 0 sinx clx Use SIMPSON s rule with 5 function evaluations to obtain an estimate S4 of the integral 0 sinx clx RICHARDSON s extrapolation suggests that B4 0694 7 S2 may be a better estimate than S4 Compute B4 and the error in each of S4 and B4 Problem 188 A Use SIMPSON s rule with 2 subintervals 3 function evaluations to es timate 02 954 also B Use SIMPSON s rule with 4 subintervals to estimate the same integral C RICHARDSON s extrapolation suggests that B4 16S4 7 S2 may be a better estimate Com pute B4 D Compute the actual error in each of the above parts Problem 189 Denote by Tn the compound trapezoidal rule with n subintervals for some function f on some interval 1 b If n is even we can apply RICHARDSON extrapolation to obtain SIMPSON s rule Sn ng 7 ngQ If n is divisible by 4 we can extrapolate once again to obtain BOOLE s rule Bn 5 7 Sn2l Given T1 1 1063582750 T2 0 8005457781 T3 0 7517873921 T4 0 7354090855 T5 0 7279356090 compute B4 2 Problem 1810 For a certain function f we have trapzoidal estimates of f dz as 1 follows T2 0 33189 T4 035947 and T3 036648 Use extrapolation to compute the SIMPSON s estimate S4 Problem 1811 Let Tn f denote the trapezoidal quadrature formula for the function f with n equal length subintervals In class we saw if n is even and we perform a certain Richardson extrapolation we obtain SIMPSON s rule in the form 4 1 Sm Tnf 7 7 420 3 3 If n is divisible by 4 we can repeat the process and we obtain BOOLE s rule in the form 16 1 Bnf Esnlf ESnQUC 1 Suppose x dx 0 4997992691593999278258 and supposeT1f 0 85826285T2f 0 057944883 and T4f 051924378 Part A Compute S4f and B4 Part b Compute and compare the errors in T4 S4 and B4 b Problem 1812 If we approximate the integral x dye of a certain function f on the interval 1 b by the compound trapezoidal rule with n subintervals T7 for n 12 and n 24 we nd T12 356782 T24 382993 Use this information to compute the SIMPSON s rule estimate Sn with n 24 subintervals Ans 348045 Mm 351 33 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis 19 GAUSSLEGENDRE Quadrature Problem 191 GAUSS quadrature on 71 1 With three nodes is given by Gem 3f flt0gt f Use G3 to estimate 2 x6 9 793 238 462 643 383 279 502 884 197 955 de Ans 2 Problem 192 The 37point GAUSS quadrature method normalized for the interval 0 4 is given by Gsf190flt2e gt f2193flt2 m Without doing too much calculation compute G3 955 bia Problem 193 Let f be continuous on 1 b n 2 1 and h k 0 1 n The compound trapezoidal rule is given by LetacC akhfor Tm g new 2fx12fx2 2fz1fzgtl Simpson s rule n even is given by 1 Snf 4TnfTgf 3 2 mo 4fx1 2M2 41 wheel he Part A Use SIMPSON s rule With 2 subintervals to estimate 02 x5 clx Part B Use SIMPSON s rule With 4 subintervals to estimate the same integral Part C RICHARDSON s extrapolation suggests that BODE s rule B4 0694 7 S2 may be a better estimate than S4 Compute B4 Part D Compute the actual error in each of the above parts and comment on it Part E Using only mental calculation deduce the approximation provided by GAUSS quadrature With 3 nodes Explain Problem 194 GAUSS quadrature With 3 nodes rescaled to the interval 0 1 is given by 5 57m 4 1 5 5 Cf flt 10 gt flt2gt Eflt 10 Find Gz5 953 without doing much computing think rst Simplify Ans 20 Miscellaneous Quadrature Problems Problem 201 We know 3 3 5 76 m for 9520 Mm 351 34 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Use this inequality and only this inequality to obtain an estimate of 1 sinx dx 0 x With an error no bigger than 10 3 Be sure to explain What you are doing and to back up your claim that the error is no bigger than 10 3 Problem 202 Determine a b and c such that the quadrature formula 00 af0bf cf1 is exact for the integral 1 x3 dx if x is a polynomial of degree 3 2 that ios1is exact for x 1 x x and x x2 Find the errorin Qx4 as an approximation to x7 dx 0 Problem 203 Determine c1 CQ Cg and C4 such that the qudrature formula QU Cif062f1csf3C4f4 is exact for the integral 00 e71 dx 0 if x is a polynomial of degree 3 3 that is is exact for x 1 x x x x2 and x3 Problem 204 Consider the quadrature method 7 5 1 fi 4 1 5 1 10 WET 0 1 1 for approximating the integral f dx Part A Given the table 0 4 1 5i 3 1 1 iim l i iily 2 E 2 10 2 100 25 2 1 5i 3 2 1 5 iff 7iim i i Einm 2 V5 5 10 2 E 40 1000 3 6 M i i i gt 63m V5 2 E 250 1000 compute the exact error in I xp for p 0 1 2 3 4 5 6 Simplify your calculations by noting some cancellations Part B If Px is a polynomial of degree 3 m Whatis the largest value of m for Which you can guarantee that P Px dx thatis for Which P is exact 0 Mm 351 35 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis 21 Iterative methods for xed points for nonlinear systems Si multaneous and successive displacements If we have a xedpoint problem for a system of nonlinear equations for example 961 91961962963 962 061962963 963 061962963 then we can try a xedpoint iteration with simultaneous displacements It glltx quotgtx gtz gtgt x9 92xl x quotaw I 93z gtz gtz gt or a xedpoint iteration with successive displacements xi 91xl x ax quot 965quot 92x 1x quotx 96 ya 96 965quot avg It should be clear how these ideas extend to systems of more equations Thprp an Amp nmmmmte afor analnnnn tom w hypothesis in the scalar case The ideas described above are actually applied to linear systems even though there are deterministic ways to solve linear systems The reason is that for very large sparse systems the iterative methods re main practical whereas the deterministic methods become unmanageable For the sake of description let s use a small system to illustrate The system of linear equations 11196 b1y 612 61 1295 2 ng 52 1396 bay 632 63 may be converted to a xed point problem by pulling out the diagonal terms 1 i b 7 x a1 61 13 6196 l y I 52 7 1295 7 92 2 l 2 753 7113957 bgyl Cs Mth 351 36 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis If we apply the method of simultaneous displacements to this xedpoint problem then we obtain the GAUSSJACOBI iteration for the original linear system If we apply the method of successive displacements then we obtain the GAUSSSEIDEL iteration for the linear system Note there is an other way to describe the GAUSSSEIDEL iteration in terms of triangualr matrices which is more convenient for theoretical work but is equivalent to our description here Problem 211 The xed point problem xfxy y 906 where x2 3y32y24y1 fx 3 39572952 4y72y23 995 y f has the solution E y 0 1645436899 02966770274 Thereis another solution 72337069 72441771 but it does not concern us here Use the initial point 950 yo 0 0 and compute 8 xed point iterates yon yn by the method of simultaneous displacements At each step compute also the error 1 2 E 7 9572 g 7 yn2gt Do you appear to have convergence Problem 212 Repeat the previous problem but this time use the method of successive dis placements Compare the errors to the ones observed in the previous problem Problem 213 Consider the xedpoint problem 3x376x224xxy753 6x2712x71y 33952 7 24y3x2 247 3095 6x2712x71y Start at the initial point 72 1 and compute 3 iterations using the method of simultaneous displace ments Your iterates should converge towards 72 92173131345952223 1 60954160414971480 Compute the 2norm of the error in each iterate Use a spreadsheet or other appropriate tool Problem 214 Redo the previous problem but this time use the method of successive dis placements Comments Problem 215 The system of equations as x27xyy2 y x2xyy2 has two solutions 0 0 and 02733011741 06027847152 Experiment with the xed point iter ation to see what happens Convert the equations into a root nding problem and experiment with NEWTON s method This problem is a bit vague Mm 351 37 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Problem 216 The system of equations 12 l 9 7 3 3y 3 1 4 1 y 77x27y277y77 has two solutions One of them is 570592062563225385 7 193100442914760222 Use the initial point 06 702 and the implied xed point iteration with the method of simultaneous displacements with 5 iterations to approximate a solution Compute the 2norm of the actual error in each iteration Comment Problem 217 The system of equations 1221 95739 333 7 121241 y Exl y y i has two solutions One of them is 570592062563225385 7 193100442914760222 Use the initial point 06 702 and the implied xed point iteration with the method of successive displacements with 5 iterations to approximate a solution Compute the 2norm of the actual error in each iteration Comment Compare the errors here with the errors in the previous problem 22 Vectors Inner Products and Norms Problem 221 Assume the canonical inner product in R5 and compute the angle in radians between the vectors 127313 and 145513 Problem 222 A student playing around with some vectors 3k k 1 2 3 4 was amused to nd 1391 2132 3133 454 4131 3132 2133 54 Are the vectors 1391 1392 1393 1394 linearly independent Justify your answer carefully in English Problem 223 Let 1 OHMDJ W1 l HWM 1 H Urqu Mm 351 38 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Part A Write 774 as a linear combination of 771 772 773 explicitly Part B Show that 771 772 773 are linearly independent roblem 224 Let 771 777 be vectors in R Let 77 be a vector in R and suppose w span771 77n The symbol means is not an element of Part A If 0777 1771 n77n 0 give a careful explanation of why it follows that a 0 Part B Now show if the set 771 777 is linearly independent then the set 771 777 77 is linearly independent 2 71 2 3 Problem 225 Consider the vectors 771 3 772 74 773 3 and77 2 For 4 1 a 9 what values of a is 77in the span of 771 772 773 2 71 2 71 Problem 226 Consider the vectors 771 3 772 74 73 3 and77 0 4 1 a 1 For what values of a is 77in the span of 771 772 773 Problem 227 There are several ways of establishing that any 5 vectors in R4 are linearly dependent Provide one or more arguments 23 Matrices Problem 231 Let C be the matrix 3 71 c l 1 1 l Find a monic polynomial p of degree 2 such that pC 0 that is nd constants a and I such that 02 a0 M 0 Problem 232 Let C be the matrix Note 25 22 7 42 423 422 122 321652 7 8 Use this fact to compute C5 without computing any powers of C You will need to do the previous problem rst Problem 233 If A alj is an n X n matrix the trace or spur trA of A is de ned by mi Z aquot 1 2 30 A4 75 and B21i Mth 351 39 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Compute AB BA trAB and trBA Do you care to make a conjecture Problem 234 Let A be a 2 X 2 matrix such that AliH 43111 1 Compute the product Problem 235 Let A be a 2 X 2 matrix such that 0 3 1 2 AMHI and411F121 Problem 236 Let A be a 2 X 2 matrix such that 40131 ABLE Find A Compute the product ts r v NH r DJ Problem 237 Let A be a 2 X 2 matrix such that 1 3 A 1 1 1 1 2 Problem 238 Let A be an n X n matrix Then the series Find A 00 tk tA 7 1c e 7 EA k0 Where A0 I converges Let 0 1 1 A 0 0 1 r 0 0 0 Compute e A Problem 239 Let 1 0 0 71 I0 I and J1 0 Mm 351 40 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Let a a and b be scalars Compute and express in terms of I and J the matrix aI bJaI 51 Problem 2310 LetA l1 2 2 3 andletBbea2 gtlt 2matrix If ABBI Where I is the 2 X 2 identity matrix compute B Problem2311 LetA 1 J First compute A3 2 A Then compute A A3 2A2 17I1423346679 Where I is the 2 X 2 identity matrix 24 GAUssJORDAN Row Reduction Problem 241 Find the complete row reduced row echelon GAUSSJORDAN canonical form of the matrix A henceforth abbreviated rref A where D5 gt 1 A5 0 Ougtm 1 0 0 Problem 242 If UCJ H ooowo oomom Or ODJO HOMOO Or OOO then nd rref B Problem 243 Consider a system of linear equations A9 After an immense number of elementary row operations we reduce the augmented matrix A b to the row reduced echelon form 0 1 0 0 3 2 0 3 0 0 1 0 2 1 0 1 nefltAEgt 0 0 0112 0 2 0 0 0 0 0 0 1 3 0 0 0 0 0 0 0 0 Call the variables in the system of linear equations x1 x2 Part A Which columns are pivotal Part B Which columns are free Part C Write the general solution of the system in vector para metric form Mth 351 41 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Let A be an m X n matrix Obviously elementary roW operations do not change roWA though they may change colA However since dimroWA dim colA we see that elementary roW operations do not change dim colA Problem 244 Consider the vectors 8 2 1 71 a 78 a 74 a 2 a 14 1 7 21 a 2 3 a 3 6 a v4 12 711 1 77 723 Find the dimension of the subspace oflR4 spanned by these vectors Problem 245 Part A Find the roW reduced echelon canonical form rref A of the matrix 2 1 7 5 4 A 1 1 5 3 3 72 2 2 72 3 Part B Identify the pivotal columns Part C What is the rank of A Part D Are the roW of A linearly independent Part E Are the columns of A linearly independent 25 LU Facrorization Problem 251 A certain 3 X 3 matrix A has an LU decomposition With 1 0 0 3 2 5 L 1 10 U 0 72 1 g 7 1 0 0 75 Solve the system Ax bWhere b 5 15 10T Answer 9 76 72 Problem 252 The 4 X 4 Hilbert matrix A given by 1 1 1 1 4 row on p m has an LU factorization A LU Where 111 1 1 11 O m Mm 351 42 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Compute L recall it is normalized 26 Positive De nite Matrices The n X n matrix A is said to be positive de nite if it is symmetric and if xTAx gt 0 for each nvector x y 0 One can show that an n X n symmetric matrix A is positive de nite if and only its principal cofactors are positive Problem 261 Show in two different ways that the matrix 2 71 A l 71 2 l is positive de nite Problem 262 Find a positive de nite matrix B such that B2 A where Alj 27 GAUssJACOBI and GAUSSSEIDEL iterative methods for lin ear systems For systems of linear equations A9 gwhere A is a sparse matrix iterative methods are often used If A is a very large matrix storage space constraints may force the use of iterative methods Consider the system of linear equations A9 gwhere the coef cient matrix A is a n X n matrix Iterative methods generally consist of decomposing the matrix A as A Q R where Q is nonsin gular and is the coef cient matrix of a system of linear equations which may be ef ciently solved for example a triangular or a diagonal matrix The corresponding xed point iteration is then Qgng 7 R550 If we take Q to be the diagonal part of A we obtain the GAUSSJACOBI iteration which converges if A is strictly diagonally row dominant If we take Q to be the lower triangular part of A we obtain the GAUssisEIDEL iteration which converges for example if A is diagonally dominant and positive de nite with positive diagonal elements The GAUssisEIDEL iteration is equivalent to doing the GAUSSJACOBI iteration of A as described above except we use the updated value of each component as it becomes available rather than Mth 351 43 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis waiting until all the new components of the new iterate have been calculated Thus the GAUssi SEIDEL method is just the method of successive displacements applied to the xed point problem obtained by pulling out the diagonal The GAUSSJACOBI method is just the method ofsimultaneous displacements applied to the xed point problem obtained by pulling out the diagonal Intuition suggests that GAUssiSEIDEL should be the better method but that is not always the case Explicitly the naming convention appears to be 1 xed point problems for systems linear or nonlinear a method of simultaneous displacements JACOBI iteration b method of successive displacements SEIDEL iteration 2 systems of linear equations a GAUSSJACOBI 2 convert to xed point problem as above and apply JACOBI iteration b GAUssiSEIDEL 2 convert to xed point problem as above and apply SEIDEL iteration Problem 271 Consider the linear system A9 Ewhere 74 71 2 a 2 A 2 5 2 and b 4 2 2 75 3 Starting with an initial guess 0 0 OF compute the rst 5 GAUSSJACOBI iterates Compute also the exact solution any way you like and use it to nd the 2norm of the error in each iterate Do the GAUSSJACOBI iterates appear to be converging As a sanity check on your work the 2norm of the error in the third iterate appears to be about 006180215 Problem 272 Consider the linear system A9 I where 74 71 2 2 A 2 5 2 and b 4 2 2 75 3 Starting with an initial guess 0 0 OF compute the rst 5 GAUssiSEIDEL iterates Compute also the exact solution any way you like and use it to nd the 2norm of the error in each iterate Do the GAUssiSEIDEL iterates appear to be converging As a sanity check on your work the 2norm of the error in the third iterate appears to be about 003885003 Problem 273 Consider the linear system A9 I where ii Hit Starting with an initial guess 0 OF compute the rst 5 GAUSSJACOBI iterates Mth 351 44 Summer 2004 B E Petersen Version 20040806 Elementary Numerical Analysis Compute also the exact solution any way you like and use it to nd the 2norm of the error in each iterate Do the GAUSSJACOBI iterates appear to be converging As a sanity check on your work the 2norm of the error in the third iterate appears to be about 018034689 Problem 274 Consider the linear system A9 I where 5 2 a 3 A2 4 and b74 Starting with an initial guess 0 OF compute the rst 5 GAUssisEIDEL iterates Compute also the exact solution any way you like and use it to nd the 2norm of the error in each iterate Do the GAUssisEIDEL iterates appear to be converging As a sanity check on your work the 2norm of the error in the third iterate appears to be about 00116276 Problem 275 Consider the system of linear equations on 2952 2953 2 73951 2952 x3 3 2952 2953 1 By separating out the diagonal convert this system to a vector xed point problem Perform one GAUSSJACOBI iteration starting with an initial guess we 0 0T What is your result Ans 13112 Problem 276 Consider the system of linear equations on 2952 2953 2 73951 2952 x3 3 2952 2953 1 By separating out the diagonal convert this system to a vector xed point problem Perform one GAUSSSEIDEL iteration starting with an initial guess we 0 0 0T What is your result Ans 192 712 Copyright 2004 Bent E Petersen Permission is granted to duplicate this document for nonipro t educational purposes provided that no alterations are made and pro vided that this copyright notice is preserved on all copies Bent E Petersen 24 hour phone numbers Department of Mathematics of ce 541 7375163 Oregon State University Corvallis OR 973314605 fax 541 7370517 email petersen mathoregonstateedu web httporegonstateedqueterseb Mth 351 45 Summer 2004

### BOOM! Enjoy Your Free Notes!

We've added these Notes to your profile, click here to view them now.

### You're already Subscribed!

Looks like you've already subscribed to StudySoup, you won't need to purchase another subscription to get this material. To access this material simply click 'View Full Document'

## Why people love StudySoup

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "I bought an awesome study guide, which helped me get an A in my Math 34B class this quarter!"

#### "I was shooting for a perfect 4.0 GPA this semester. Having StudySoup as a study aid was critical to helping me achieve my goal...and I nailed it!"

#### "Their 'Elite Notetakers' are making over $1,200/month in sales by creating high quality content that helps their classmates in a time of need."

### Refund Policy

#### STUDYSOUP CANCELLATION POLICY

All subscriptions to StudySoup are paid in full at the time of subscribing. To change your credit card information or to cancel your subscription, go to "Edit Settings". All credit card information will be available there. If you should decide to cancel your subscription, it will continue to be valid until the next payment period, as all payments for the current period were made in advance. For special circumstances, please email support@studysoup.com

#### STUDYSOUP REFUND POLICY

StudySoup has more than 1 million course-specific study resources to help students study smarter. If you’re having trouble finding what you’re looking for, our customer support team can help you find what you need! Feel free to contact them here: support@studysoup.com

Recurring Subscriptions: If you have canceled your recurring subscription on the day of renewal and have not downloaded any documents, you may request a refund by submitting an email to support@studysoup.com

Satisfaction Guarantee: If you’re not satisfied with your subscription, you can contact us for further help. Contact must be made within 3 business days of your subscription purchase and your refund request will be subject for review.

Please Note: Refunds can never be provided more than 30 days after the initial purchase date regardless of your activity on the site.