### Create a StudySoup account

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

Already have a StudySoup account? Login here

# COMPUTATIONAL METHODS MCEN 3030

GPA 3.73

### View Full Document

## 146

## 0

## Popular in Course

## Popular in Mechanical Engineering

This 30 page Class Notes was uploaded by Deron Schmitt on Thursday October 29, 2015. The Class Notes belongs to MCEN 3030 at University of Colorado at Boulder taught by Oleg Vasilyev in Fall. Since its upload, it has received 146 views. For similar materials see /class/232013/mcen-3030-university-of-colorado-at-boulder in Mechanical Engineering at University of Colorado at Boulder.

## Reviews for COMPUTATIONAL METHODS

### 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/29/15

Chapter 3 Interpolation Interpolation is de ned as a continuous representation of a discrete set of data points Applications might be for differentiation or integration or simply estimating the value of the function between two data points There are two distinct interpolation methods to obtain such a representation 0 standard where the curve passes smoothly through the data points and 0 least squares where a smooth curve passes suf ciently close to data points that have some uncer tainty Note that with standard interpolation we actually pass a curve through the data If data is from an experiment with some uncertainty it is best to use the method of least squares which does not require the approximating function to pass through all the data points Our use of interpolation will be primarily in deriving differentiation and integration formulas thus we will use standard interpolation exclusively When the interpolating function is evaluated outside the de ning range then it is called an extrapolation 31 General Interpolation Problem The general standard interpolation problem is de ned as follows Given 1 a set of data points zii 01 N 2 a set ofvalues E 01N and 3 a set of basis functions bjzj 01 N Find coef cients ajj 01 N so that N Zajbjtzi 239 01 N j0 Note that for simplicity we have assumed linearity in the coef cients aj and subsequently the above equation is just a system of linear algebraic equations which can be written as f i Bijajv where Bij E 12111 The set of basis functions bj is still arbitrary at this point as long as they are complete in the space of f but for convenience we impose two restrictions on them 1 they should be easy to evaluate and 2 they should be easy to differentiate and integrate The most common basis functions used that satisfy these restrictions are polynomials of various kinds 13 14 CHAPTER 3 IN TERPOLATION 32 Lagrange Interpolation Suppose we have a set of N l nonequally spaced data By choosing the basis functions bj zjj 01 l i i N we can construct a polynomial of degree N that passes through the data N V 2 N 125 ajzg a0alziagzi aNz 20lNi 10 Here we have N 1 equations for the N l unknowns a0 a1 aNi This procedure for nding the coef cients of the polynomial is not very attractive It involves solving a system of algebraic equations which is generally illconditioned for large N the resulting matrix is called a Vandermonde matrix In practice we will de ne the polynomial in an explicit way as opposed to solving a system of equations De ne a polynomial of degree N associated with each point zj L Ngtltzgt aw i zogtltz e z1gtltz 7 some 7 1H1 z 7m where aj is a constant to be determined Let7s use the product notation for LEN N N L z aj H I 7 20 2 If I is equal to any of the data points except zj LEN 0 Note N 119011 11 H 11 Ii 1 0 2 Let 1 N 11 HW M 7 3 so that LNI I i 10 quotI i 11701 1H1 39 39 39 I IN J Ijirol39lrjirjemrj Ij1quot39rjIN Then N L n 0 ifia jy L1 Il 6u 1ifijm where llj is the Kronecker delta function Next form the linear combination N fltzgt 4 sz This is a polynomial of degree N because it is a linear combination of polynomials of degree Ni It is the classic Lagrange polynomial By construction it passes through all the data points For example at z zl fltzigt L3Ngtltzigtfo L Ngtltzigtf1 L Ngtltzigtfj L gtltzigtfNi Since LEN is equal to zero except for j i and LENzZ l flquot Thus is the desired polynomial Note that the polynomial interpolation is unique That is there is only one polynomial of degree N passing through a set of N 1 points The Lagrange polynomial is 33 PIECEWISE LAGRAN GE IN TERPOLATION 15 just a compact and numerically better behaved way of expressing the polynomial whose coefficients could have also been obtained from solving a system of algebraic equations If we assume the data are exactly the values a smooth function takes at 11 z39e then it can be shown that MN 1 5 N 1 Note that Fz z 7 10 I 7 IN and 07 so that the error is largest between points and for more wiggly functions For a large set of data points say greater than 107 polynomial interpolation can be very dangerous Although the polynomial is fixed tied down at the data points7 it can wander wildly between the points which can lead to large errors for derivatives or interpolated valuesi Better results are obtained by piecewise Lagrange interpolationi Example Take e for 0 S I S 1 We use three points N 2 IO 07 11 0512 17 and take Then WE fr FI7 IOlt5ltINA 7 I MHI W I 1 10II2 I 1 10I11 I flag 7 10 7 11 10 7 I2 0 T 11 7 mo 11 7 I2 1 T 12 7 10 I2 7 11 2 At 1 025 we have 70 2570 75 0 025 7075 05 025 70 s 71 W H 05 705 W 0 37500012365417 13397857 1271756 025 7025 1 lt1gtlt05gt The exact value of the function at z 025 is y0i25 603925 12840254 so that the error is approxi mately 1i 33 Piecewise Lagrange Interpolation Instead of fitting a single polynomial of order N to all the data7 one fits lower order polynomials to sections of it This is used in many practical applications and is the basis for some numerical methods The simplest case is linear Lagrange interpolation z 1 7 z z 7 z 3 71 fi 71 fir n1 Ii n1 Ii The subscript on indicates that this formula holds only for 11 S I lt 1141 An improvement can be obtained by using the quadratic Lagrange polynomial I I N I 7 ziz 7 n1 f1 1171 7 ziri71 EH1 z 7 Ii71z Ii n1 7 Ii711i1 7 xi fi71 7 Iiilxz 7 n1 f39 z 7 Ii711i 7 n1 Z fi17 for 1171 S I S 1141 This procedure can be extended to higher order interpolationi 34 Cubic Splines The main problem with piecewise Lagrange interpolation is that it has discontinuous derivatives at the boundaries ofthe segments which causes dif culties when differentiating lnterpolation with cubic splines1 circumvents this dif culty e define a spline interpolating curve by the following criteria 1The equation from mechanics describing the shape of the drawing tool called a spline is EIfi m F where E is the modulus of elasticity I the moment of inertia and F the applied axial force In our case 0V 16 CHAPTER 3 IN TERPOLATION H The curve is piecewise cubic that is the coefficients of the polynomial are different on each interval Ii1i1 i to i The curve passes through the given data 9 With the exception of the end points the first and second derivatives are continuous at the points Iii From continuity of f z using the piecewise linear Lagrange formula we have zi1 I I Ii l E l 39 Jr 1 39 i L i 1 i l i 1 i 11 Integrating the equation twice and tting to the data and fzi1 results in mamgt3 i 7 3 hi mac f ltz1gt 6h flt1i11i1 1 h gf xi fin 7 f zi1 7 where hi 14171 If we can find f the formula will be complete This is obtained from continuity of f using two adjacent formulasi The resut is hi4 hi f1i1 7 i 1121 3 hi i h1 h LINIiil fIi if zi1 6 6 13971 This is a tridiagonal system for f zi that is numerically easily solvedi This equation cannot be applied at i 0 or N so that we only have N 7 1 equations for N 1 unknowns We get two additional equations by making some assumptions andor approximations in the end intervalsi Some choices of boundary conditions are 1 Periodic MID f zN1fIN my 2 Paralmlz39c MID my MIN MINA 3 Natural or free MID MIN 01 4 Cantileven MID Am MIN Af IN1 0 g A g 1 For equally spaced intervals h we can show that if are the exact values a smooth function takes at Ii t e h4 m 7 ms 4316 re lt 5 lt INT Example The simplest case requires 4 points N 3 two in the middle and two end points Also we assume equally spaced intervals and cantilever end conditionsi Then for i 1 Wm 4f I1 f I2 Waco 7 2f11 fr2l 6 64fl 10 f0 2161 f2 Applying the end condition at i 0 6 4 Mfl f f0 2fl f2 Similarly at i 2 6 i 4Af f1 2f2 fs 35 PARAMETRIC IN TERPOLATION 17 These equations can be readily solved for ff and To make things simple we will evaluate the function in the middle of the interval at z 11 zg2 so h2 7 ltf1f2gt7gltfrfrgt7 3 f1 f2 7060 fs lisltsiiAgtl 8lt5Agt For an 91139 5117 0 S I S11 A0 f0i57e 5 m5 X10 3 1 f0i5 7 eO m 5 gtlt10 4i In most cases the spline gives a smooth t to the data In some cases however it wiggles too much A cure for this is to put the spline in tension The effect of tension is to replace equation 31 by 1 mac 7 02121 7 Wm 7 WW 7 I1 1 139 i1 Ii lf 1i1 02f1i1l7 On integrating twice and requiring that the curve pass through the points we obtain 7 sinh lama Il f 1i n1 1 Wm 1611 7 sinh azi1 7 a2 T Ii1 7 7 02 T sinh 01 7 f zi1 I 7 sinh azi1 7 a2 T I 7 z fzi17 139 Then on applying continuity of the rst derivative we nd the equation for f to be 1 a f zi1 acoshahFl 1 acoshahi 7 i f zi 02 sinh Hill1 hpl sinh ahi hi 02 i 7 a f 1i1 f1i1 7 10 f1i71 hi sinhahi hFl sinh Hill1 02 z 13971 Again we have a tridiagonal system of equations which is easily solvedi Generally it is recommended that if the natural spline is not satisfactory one should try the tension spline with the tension parameter set equal to unity as a rst guess If the results still wiggle too much the tension 0 can be increased Endpoint conditions for the tension spline may be modeled after those for the natural spline by replacing f zi by f zi 7 02fzi in themi 35 Parametric Interpolation None of the methods discussed is capable of interpolating curves having in nite slopes such as in com puting contours fz y constanti One way out of this dif culty is to use parametric interpolationi We regard both I and y as dependent on a third variable 8 The choice of the parameter s is free as long as it increases monotonically along the curve A simple choice is the arclengthi Then we take 8 so 0 S S1 80 11 i 10V 91 yo2l127 S S2 81 12 7 102 y2 7 MVP27 and SO 0117 The parametric curves 18 and ys are treated by the usual spline method For a closed curve periodic end conditions are appropriate for other cases other end conditions may be applied 36 Multidimensional Interpolation The simplest methods in multidimensional interpolation are extensions of onedimensional methods For example if a function of z and y is given on a regular array of points then one can do the following For each line I I on which the data are given compute the second derivatives which we call of 18 CHAPTER 3 IN TERPOLATION the spline t in y by the usual method The second derivatives then may be considered functions of z and may be interpolated by the spline tting in 1 To evaluate the interpolated value at a given 17y7 we rst determine in which box it is placed Given that 11 lt z lt 1 and yj lt y lt yj17 we can apply the spline ts to f zi7 yj and f zi1yj to nd f zyji In the same manner7 from f ziyj1 and f zi1yj1 we can nd f z7 yj1i Finally7 using f zyj and f zyj1 we can compute Lectures 5 t0 7 Direct methods for solving Ax b Notation A E Rmxm and Lb E R a11 a12 quot39 quot39 alm b1 1 a21 a22 02m b2 x2 A l 133 7 b 7 i bmil zmil am amz amm bm zm Ax b represents a system of linear equations to solve for x a11 112 aim m b1 a21 122 12me 52 am i am z amm m bm Gaussian elimination In the above system it is possible to eliminate all the aklzl by multiplying the rst line by 73 and subtracting it to the k th line Doing so will lead to the modi ed system a11 012 aim m b1 0 192 1me by 0 05332962 alianm big The modi ed values 11 and b0 are resulting from the subtraction an i a a an ij i ii Iii a11 a 1 bk 7 blai 2 where 2397j gt 1 and k gt 1 The step of eliminating the am for k gt 1 in the system can be written as a matrix operation called a Gaussian transformation 1 0 0 0 fanan 1 0 0 iaglan 0 1 M1 iamlan 0 1 or equivalently as M1 I 1 777115 where I is the m gtlt m identity matrix7 elT 1000 70 and 771 07702111177 1310117 7am1011l us7 0 a atquot 0 a1 a1 MA 22 2 lo ax aw 51 by M11 b9 b1 The second elimination step is then written as MZMlAz MZMlb 4 where the Gauss elimination matrix is 1 0 00 0 1 0 0 M2 0 fagall 1 16 milas o 1 11 and a11 a12 alm 0 0 0 MngA 3 0 1g ang a a gig 032 The product of the 2 transformations is 1 0 0 0 fanan 1 0 0 M2M1 031 111 faga9 1 7am1a11 fagga lz 0 1 In general the Gauss transformation is noted as where 6k is a vector of zeros with a one at the k th place and 0 0 Wk kil kil aiwikvaigk kill kil aink Valsk is the Gauss vector having the rst k places set to zero They also enjoy the nice property Mle mlelTI mkeg I mlelT mkeg mleleke T T I mlel mkek if and only if l S k This is due to the fact that mk has its rst k elernents set to zero Thus7 elek will be zero only if the only one in 61 is in the rst k th positions 1 S l S Carrying the procedure until the last block is reached brings to the upper triangular system Ui Mmileig 39 39 39 M2M1A Mmileig 39 39 39 MZMlb 3 or simply Um Explicitly the matrix U is a11 a12 alm 0 all a ii U 0 am a 0 0 a333 It is easy to see that zm blif la 1 since the vector is is b1 b9 6 l J In a similar fashion7 since we now know am it is possible to obtain zm1 MK 7 alime mil W72 amil nil The multiplications of A by the matrices Mk transforming A to U is the Gauss elimination step Solving for z with the matrix U is the backwardsubstituti0n step It is possible to write a mathematical formula for the backsubstitution bm kk l 7 7quot k altmik7mw m7 m7 1 mik 7 957M W 7 amikynik The pseudo code performing the transformation of A into U is given a b and 111 do k1m 1 do ik1m fac aikakk do jk1m aij aij fac akj end do bi bi fac bk end do end do the output ofthe above code returns U inside A the upper triangular part of A is replaced by U and b The pseudo code performing the backsubstitution is simply do kmm 11 xk bk do ik1m Xk Xk aki Xi end do xk Xkakk end do Up until now we have assumed that the divisions by diagonal elements during the elimination process was de ned all non zeros In the case of a zero diagonal element one needs to pivot lines or rows in the matrim The subject will be discussed later LU factorization Gaussian elimination is not very attractive when the system needs to be solved for multiple right hand sides rhs This is obvious since b is trans formed to l using the Gaussian transformation matrices Mk This change of b is also depicted in the Gaussian elimination pseudo code In what follows we will show that in fact the inverse of the combined Mk putting A in upper triangular form is a lower triangular matrix that we will call L Generating L requires only a simple modi cation to the Gaussian elimination pseudo code Thus the goal of this section is to transform A into LU or in other words construct the identity A LU This is quite useful if the following steps are considered 0 Factorize A into LU 0 Given I solve Ax LUx b 7 Solve Ly b using forward substitution 7 Solve Us y using backsubsttutton The last two sub steps are repeatable for different Us and can reuse the factorization of A into LU Mathematically and straight from the Gaussian elimination we have ob tained Mm1Mm2 M2M1A U inverting the Gaussian transformation matrices leads to A MflM1M 2M 1U Now if it is possible to show that Mf1M 1M 2M 1 is a lower triangu lar matrix we are done To do so lets multiply one Gaussian transformation with a tentative inverse I 7 mke Mk I 7 mke I mkeg I 7 mkeg mkeg 7 mkeg nkeg I 7 mkeg nkeg I Notice that the last piece is zero since 6k has one at position k and zero elsewhere while mk has its rst k positions 123 k set to zero their multiplication egmk gives 0 E R We have thus shown that the inverse of a Gaussian transformation is a Gaussian transformation with a change of sign The lower triangular matrix can be constructed noticing that m71 71 71 71 71 i T M1 M2 Mm72Mm717I7 E mlel 11 this property was shown in the previous section for the Gaussian elimina tion In matrix form we obtain 1 0 0 0 021011 1 0 0 MflM1M12M 1 031011 ag12a212 1 39 laman gasas 7aii a ll The Gauss elimination pseudo code can be modi ed to create the L and U matrices Without any pivoting strategy7 we are going to call this algorithm naive LU given a andm set L 0 do k1m 1 do ik1 m fac ai1akk Lik fac aik 0 do jk1m aij aij fac akj end do end do end do do k1m Lkk1 end do at completion7 A is transformed into U with zeros in its lower part and L is created Pivoting strategies The pivot is the element of the matrix with whom we perform divisions in both the Gaussian elimination and the LU factorization A zero pivot obvi ously pose problem All hopes are not lost since if A is regular then a non zero pivot must exist in the same column 1 Partial pivoting before eliminating values in column k nd the largest lagginl for 239 k k 17m The largest value becomes the pivot instead of i k4 the naive akk 2 Full pivoting before eliminating values in column k nd the largest maxla 71lla 1l for 239 kk 17 m The largest value becomes the i i i kil pivot instead of the naive akk We will interest ourselves to the rst approach only since it involves only row interchanges while the second performs both row and column inter changes LU With partial pivoting given a and m set L 0 do k1m 1 q findpivotakk tmp Pk Pk Pq Pq tmp Exchange row k with row q in a Exchange row k with row q in L do ik1 m fac aikakk Lik fac aik 0 do jk1 m aij aij fac akj end do end do end do do k1 m Lkk1 end do The pseudo code performing the forward and backward substitution is simply given P L U b and 111 do k1 m bnewkb POO end do do k12 m yk bnewk do i1k 1 yk ykLki yi end do end do do kmm 1 1 xk yk do ik1 m Xk XkUki xi end do xk xkUkk end do The resulting factorization is PA LU implying that A PTLU The dif cult part is to modify L Consider the matrices Pk representing the permutations of the rows identity matrix with the column interchanged happening at each step of the factorization then Mm1Pm1M2P2M1P1A U 8 The inverse of the above leads to A PlTMflPZTMZ l Pg1M 1U PA Pm1Pm2P2P1A PPlTMflPZTMZ l Pg1M 1U by inspection we can see that L I Pm2Pmg P1Mf1 7 I Pm2Pmg P2M2 1 7 I Pm72PmesPsM 1 I M111 1 9 which is only the Gaussian transformation inverse with the permuta tions namely My I 7 RH Pk1Pkmke 10 with My My no permutations on the last transformation 10 Examples 1 Gaussian elimination with partial pivoting 3 17 10 A 2 4 i2 6 18 712 the pivot in the rst column is 6 0 0 1 P1 0 1 0 1 0 0 6 18 712 HA 2 4 i2 3 17 10 the Gaussian transformation is 1 0 0 M1 713 1 0 712 0 1 we get 6 18 712 M1P1A 0 72 2 0 8 16 The second pivot is 8 in the second column 1 0 0 P2 0 0 1 0 1 0 thus 6 18 712 P2M1P1A 0 8 16 0 72 2 the second Gauss transformation is 1 0 0 M2 0 1 0 0 14 1 and nally we obtain 6 18 712 M2P2M1P1A 0 8 16 0 0 6 Note that if given a b then 13 MngMlPlb 2 LU with partial pivoting Using the same matrix but not considering the ones in the Mk the permuted inverse Gaussian transformation is minus identity 0 0 P1Mf1 i I 12 0 0 0 0 proceed as before last step no reordering 0 0 0 P1Mf171M2 171 12 0 0 13 i 4 0 adding the diagonal and multiplying by U on the right 1 0 0 6 18 712 1 0 0 0 0 1 3 17 10 12 1 0 0 8 16 0 0 1 0 1 0 2 4 2 13 714 1 0 0 6 0 1 0 1 0 0 6 18 712 WithPP2P1 1 0 0 6 18 712 0 0 1 3 17 10 12 1 0 0 8 16 1 0 0 2 4 2 13 714 1 0 0 6 0 1 0 6 18 712 multiplying A with P 1 0 0 6 18 712 6 18 712 12 1 0 0 8 16 3 17 10 13 714 1 0 0 6 2 4 i2 Chapter 5 Numerical Integration In this section we will discuss methods for evaluation of the integral of the function in the interval 17 I E Ab z m The function is de ned either analytically or on a set of discrete points a 1011 i i i IN 12 All integration formulas also called quadrature formulas are of the form N I W 2115me j0 where wj are usually referred to as weights The speci c choice of these weights is what distinguishes the different quadrature formulasi 51 NewtonCotes Formulas One of the easiest ways to obtain quadrature formulas is to use Lagrange interpolation on an evenly spaced mesh and integrate the result Divide the range of integration into N intervals I 7 a N 20liuNi Pass a Lagrange interpolating polynomial of degree N through the points Ii ziaih h N Pltzgt ELENVIVW Now integrate this polynomial b N b N N I Pam me L MW b 7 a Zojmzj a F0 a F0 where b CJN bia L Ngtzdx are pure numbers called Cotes numbers and are given in the table belowi Observe that in this case the weights are wj be aCJNi 23 24 CHAPTER 5 NUMERICAL INTEGRATION It can be easily shown that N EON 1 and CJN 0 j0 N M MOO MCl MC2 MC3 MC4 MC5 MC6 Truncation Error 0 1 1 412 X10 2h f 1 2 1 1 813 X 107271310 2 6 1 4 1 315 X 10 4h5fi 3 8 1 3 3 1 116 X 107471516011 4 90 7 32 12 32 7 512 X 1077717160 5 288 19 75 50 50 75 19 316 X 1077717160 6 840 41 216 27 272 27 216 41 614 X10 10h9f i Table 511 NewtonCotes coefficients Using Table 511 we recover some well known formulas N 0 Rectangular rule 1 biafa Midpoint rule 1 b71010 2 N 1 Trapezoidal rule 1 1 bafafbl N 2 Simpson s rule 1 a b 1gbafa4flt 2 gtfltbgtl etc Usually the formulas are not applied over the whole interval but at each subinterval 239e we use piecewice Lagrange intepolating polynomials rather than a single interpolation over the complete integration interval Subsequently the error estimates are for a single interval it should be summed up over all intervals to get an overall error estimate 52 Trapezoidal Rule For one interval trapezoidal rule is given by Tz1 hi 1239 90d 306739 fi1 For the entire interval ab the trapezoidal rule is N71 1 1 N71 Igliwhlt2f02fNj2fjgt7 51 where uniform spacing is assumed 53 SIMPSON S RULE 25 53 Simpson s Rule Simpson7s rule for uniform intervals is given by h N71 N72 Iwg fofN4 Z 19 Z fj 52 71 72 PM To use Simpson7s rule for the entire interval of integration the number of points7 N l7 must be odd even number of panels i 54 Error Analysis Consider the midpoint rule for the interval zizi1 Tz1 mm 7 him where yi zi12 is the midpoint of the interval Let7s replace the integrand with its Taylor series about yi me 7 mi 7 z 7 mm 7 gm 7 yawya 7 gm 7 yawquotya so that as 1 mm 7 him 7 z 7 W2 71 mi 7 z 7 m3 71 mi 7 as All the terms with even powers of I i vanish7 and we obtain mm 7 him 7 iwwyn hwwi a 53gt Thus7 for one interval the midpoint rule is third order accurate Now let us perform an error analysis for trapezoidal rule Consider Taylor series expansions mi 7 mi 7 gillHm 7 gm 7 4178h fmltyigt zz1 7 mi 7 gillmi hgf myi 4178h fmltyigt Add and divide by two 141 7 I 1 2 1 4 iv f 7 M hif lty1gt hif lty1gt a Solve for and substitute into equation 53 TI 239 i 1 1 iv mm 7 7 Emmi 7 Wm My 7 a 54gt Thus7 for one interval trapezoidal rule is also third order accurate and its leading truncation error is twice in magnitude but of opposite sign of the truncation error of the midpoint ruler 26 CHAPTER 5 NUMERICAL INTEGRATION To obtain the integral for the entire domain we sum both sides of equation 54 assuming uniform spacing we have N71 751 h N71 hg N71 I fzdz f0 fb2fj12f yi i0 Tr j1 i0 15 N71 lt39 gt 2316 90 55 Now we will apply the mean value theorem to the summationsi The mean value theorem states that there exists a point E in the interval a 12 such that N71 2 mi we i0 Similarly there is a point Ein ab such that N71 I I Z M ya NfltwgtltsgtA i0 Noting that N b 7 ah and substituting in equation 55 we obtain h4 N71 2 1 fafb2 f1 7 biaf f 7 b afiv5 iw 5 6 Thus the trapezoidal rule for the entire interval is second order accuratei One can easily show that Simpson s formula for one panel 113 142 can be written as 5U Mltfgt gm where and denote midpoint and trapezoidal rules applied to the function Note that the midpoint of the interval is 1141 Using equations 53 and 54 and the mean value theorem we see that Simpsonls rule is fourth order accurate for the entire interval abi 55 Trapezoidal Rule with EndCorrection This rule is easily derived by simply substituting in equation 54 the second order central difference formula for f 7 fifi1 1 sfz391 fi 5 11 7 hi i 0hi Let hi h const then h N71 h2 N71 1 E 2 12 fi1gt E Z ft1 0h4 i0 i0 Cancellations in the second summation leads to h N71 hg 7 7 7 7 7 7 4 I 7 2 m 1 2 f1 12 f b f ltagti W gt Thus trapezoidal rule with endcorrection is fourth order accurate and can be readily applied provided that the derivatives of the integrand at the end points are known 56 RICHARDSON EXTRAPOLATION AND ROMBERG INTEGRATION 27 56 Richardson Extrapolation and Romberg Integration Richardson extrapolation is a powerful technique for obtaining an accurate numerical solution of a quan tity eg integral derivative etc by combining less accurate solutions The essential ingredient for application of the technique is knowledge of the form of the truncation error of the basic numerical method used We shall t t quot 39 of Richard rm A 39 by using it to improve the accuracy of the integral I Ab z z with the trapezoidal rule as the basic numerical method From our error analysis for trapezoidal rule we have g c1h262h403h6 H 57 N71 W fan 2 fj h 1175 N71 fa 7 fb 2 19 l as The trapezoidal expressions can be written as I1 7 I 7 51h 7 m4 7 53H 7 59 Now suppose we evaluate the integral with half the step size hl h2l Let7s call this estimate I2 2 4 5 I2I7c1h77c2hi h 4 16 53a 510 Now we can eliminate 0h2 terms by taking a linear combination of equations 59 and 510 4k 7 I1 7 1 4 5 5 T7IEC2h Jr cSh i 511 This is a fourth order approximation for I In fact equation 511 is a rediscovery of Simpson7s rule We have combined two estimates of I to obtain a more accurate estimate This procedure is called Richardson7s extrapolation We continue with this procedure Let7s evaluate I with h2 h4 h2 h4 hF IS1761E76225676340967m 5 12 Combining with I2 we get 4I3 7 I2 7 1 4 5 5 3 7I64C2h 1024c3h i 513 Now we can eliminate 0h4 terms between equations 5 11 and 513 The result will be sixth order accurate This process can be continued inde nitely The essence of the Romberg integration algorithm just described is illustrated in the following diagram 00 00 0h5 11 I2 Eqnl 511 I3 Eqnl 5 13 7 Note that in each step of calculations half of the needed functions have already been used in the preceding calculation and need not be recalculatedl Also in principle arbitrary accuracy can be achieved In practice we are limited by the particular machine accuracyl 28 CHAPTER 5 NUMERICAL INTEGRATION 57 Adaptive Quadrature Often it is wasteful to use the same mesh size everywhere in the interval of integration a 12 Ideally one would use a ne mesh in the regions of rapid functional variation and a coarse mesh where the integrand is varying slowlyi Adaptive quadrature techniques automatically determine panel sizes in various regions so that the computed result meets some prescribed accuracy requirement supplied by the user That is with minimum number of function evaluations we would like a numerical estimate T of the integral such that S 57 14 Abfzdzif where e is the error tolerance provided by the user but in such a way that the local error 12quot ii S 6i57 where Ii and denote the exact and approximate integrals in the interval 113 1H1 is uniformly dis tributed To demonstrate the technique we will use Simpsonls rule as the base methodi Let7s divide the interval ab into subintervals 1131141 Divide this interval into two panels and use Simpsonls rule hi hi 5 g mm W Egt m higtl Now divide the interval into four panels and obtain another estimate for the integral hi 5 5 Am 4m gt 2m gt 41 Sh 4 m my The basic idea is to compare the two approximations SE1 and SEQ and obtain an estimate for the accuracy of SE2 i If the accuracy is acceptable we will use SE2 for the interval otherwise the method further subdivides the interval From our error analysis we know that I 7 50 ch fi zi g 5 14 Sh 4 h 5 h iv Iiis 2gtclt gt flv1i1f1i Each of the terms in the bracket can be expanded in a Taylor series about the point hiZ h h h h m 101 1101 z f n4 f n2 4f zl2 Sh h h h W 71 W J 711 J f 11 4 f 1124f 112 i Thus 5 1175 25 fltivgtzi 5 15 Subtracting 515 from 514 lt2gt7lt1gt1jm E Si Si 16Chzf I 2 A 58 GAUSS QUADRATURE 29 So the error in SE2 is about 115 of the difference between 511 and 52 wai IZ SZ 15 Sf 7 Siml 1 hi b7a 1 E 5 7 S 1gtl g e 7 616 516 then SE2 is sufficiently accurate for the interval 1131141 and we move on to the next intervali 1f condition 5116 is not satis ed the interval will be subdivided further We note that condition 5116 guarantees that e is an upper bound on the total integration error since 17f 7 N71 1759 7i N71 SW7ng lt iNil SW7ng lt LN71h75 l 31 This is the essence of adaptive quadrature algorithms Similar methodology can be divised for other base methods such as the trapezoidal rule As with Richardson extrapolation the knowledge of the truncation error can be used to obtain estimates for the accuracy of the numerical solution without knowing the exact solution 5 8 Gauss Quadrature Recall that any quadrature formula can be written as b N I m a 10 1f the function f is given analytically or we are at liberty to choose the locations I we have two important choices to ma e e ave to select the location of the points zj and the weights wji The main idea behind Gauss quadrature is to make these choices for best accuracy The criterion for accuracy is the highest degree polynomial that can be integrated exactlyi Recall that the trapezoidal rule integrates a straight line exactly and Simpsonls rule integrates a cubic exactlyi In general NewtonCotes algorithms integrate polynomials up to order N N even with N 1 points exactlyi As we will show below Gauss quadrature integrates a polynomial of degree 2N 1 using only N 1 points Let f be a polynomial of degree 2N1i Suppose we represent f by an Nch order Lagrange polynomial Pi Let zozliiizN be the points on the z axis where the function f is evaluated Using Lagrange interpolation we have N Pltzgt Z L Ngtltzgtfltzjgt This representation is exact if f is a polynomial of degree Ni 7 131 is a polynomial of degree 2N1 which vanishes at 10 111 1N because P was constructed to pass through fzo fzl i i i at zozli i i 1N1 et Fz z 7 zoz 711 I 7zNi Thus we can write the difference 7 131 in the following form N 161 7 Pm 7 19024 k0 lntegrating this equation for the moment without af xing limits fzdzPzdzFzlqkzkdzi 30 CHAPTER 5 NUMERICAL INTEGRATION Interchange summation and integration7 and suppose we demand that Fzzkdz 07 k 0717M7Ni This corresponds to a system of N 1 equations for 107117 i i i 7 zNi We will choose zo7zl7i i i 7 IN such that this condition is satis ed Choosing the abscissas in this manner leads to the following expression for the integral Pzdz 210139 sz where wj LNzdz are the weights F is a polynomial of degree N 1 that is orthogonal to all polynomials of degree less than or equal to Ni zo7zl7i i i 7 IN are the zeros of this polynomiali These polynomials are called Legendre polynomials when I varies between 71 and 17 and they can be generated by Rodrigues7 formula lt71 d 2 n F7 2 m 1 1 7 where Fn is the Legendre polynomial of degree n eg7 Foz l7 F1z I7 F2z 312 7 D7 They are orthogonal to each other7 that is 1 2 F F d 76 i 1 I 2n1 77m Their zeros are documented in Table 52 Having the zeros7 the weights wj can be readily computed7 and they are also documented in Table 52 In the table N7 not N 17 denotes the number of points Note that we can always transform the interval a S I S 12 into 71 S 5 S l by the tranformation 7 b a b 7 a 17225i Typically7 to use LegendreGauss quadrature tables to evaluate the integral Ab frdx7 one rst changes the independent variable to 5 and obtains the weights wj and the points on the abscissa7 507 517 i i 7 EN7 from the tables for the chosen N The integral is then approximated by N biaZ ba bia 2 j0wjflt 2 T 2 In addition7 if one wishes to improve the accuracy by going to a higher order method7 the calculation must be done from scratch since the abscissas of the various orders are different 59 SemiIn nite Intervals If we have the integral 1 Amyrdr7 510 INFINITE IN TERVALS 31 rede ne the function so that 1 67271 thus 00 00 N N I 0 gzdz 0 e Efzdz m ijeEVg j The abscissas then turn out to be zeros of the Nch order Laguerre polynomial e d H dz 6 eg L0 l L1z 71 1 L2z 12 7 4x 2 H and have the orthonormal property Ln 11 e ELmzLnzdz Smut 0 The zeros and weights are given in Table 53 Note that in the table N not N 1 denotes the number of points 510 In nite Intervals If we have the integral 1 y a fltzgt em rede ne the function so that thus 00 oo 2 N N 2 I gzdz 6 7 m ijexig j foo 00 j0 10 The abscissas then turn out to be zeros of the Nch order Hermite polynomial 7 12 d 7052 1151 e d lt6 gt eg H0 l H1z 21 H2 4x2 7 2 i i and have the orthogonality property e EZHmzHnzdz 2nnl6mni The zeros and weights are given in Table 54 Note that in the table N not N 1 denotes the number of points 511 Singularities Two possible methods to remove singularities are presented by examples 1 Integration by parts eg 7r 7r coszdz 2 coszlg2 sinzdzi 0 f 0 I The procedure may have to repeated more than once CHAPTER 5 NUMERICAL INTEGRATION 2 Subtraction eg 7rcosxd dz 7rcosrild z i 7r o xE o xE o f Now the rst part can be integrated analytically and the second part is not singular any more 51 1 SINGULARITIES n if 1 E 1131037 at Abscissasir 1 Zeros of Legendre Polynomials 057735 000000 0774 9 033998 086113 000000 053846 090617 023861 066120 093246 000000 040584 074153 094910 iJ39i 02091 00000 66692 10435 63115 00000 93101 98459 91860 93864 95142 00000 51513 11855 9123 Weight Factorsx 7171 u if 2 728 quot 018343 46424 95650 0 36268 89626 100000 00000 00000 052553 24099 16329 031370 079666 64774 13627 022238 n3 096028 98564 97536 010122 88880 888888 88888 88889 41483 955555 55555 55556 59 000000 00000 00000 033023 quot74 032425 34234 03809 031254 84856 065214 51548 62546 861337 14327 80590 026861 94953 834785 48451 37454 083683 11073 26636 818064 096816 82395 07626 088127 It 0 88880 056888 88888 88889 7210 65683 947862 86704 99366 814887 43389 81631 029552 38664 023692 68850 56189 043339 53941 29247 026926 867948 95682 99824 021908 4 986506 33666 88985 014945 83197 046791 39345 72691 997399 65285 17172 906667 66265 936876 15739 48139 03152 017132 44923 79179 vi12 012523 34885 11469 824914 n 036783 14989 98180 023349 88888 041795 91836 73469 258731 79542 86617 020316 77397 038183 80585 05119 076990 26741 94305 016007 99394 927979 53914 89277 899411 72563 70475 810693 42759 012948 49661 68870 898156 06342 46719 004717 2176 7216 W 889501 25098 37637 448185 018945 86194 55068 496285 828160 35507 79258 913230 018268 34159 44923 588867 045881 67776 57227 386342 016915 65193 95002 538189 061787 62444 02643 748447 814959 59888 16576 732081 875548 44083 55003 833895 812462 89712 55533 872852 086563 12823 87831 743880 809515 85116 82492 784810 894457 50230 73232 576078 906225 35239 38647 892863 095940 09349 91649 932596 002715 24594 11754 094852 R20 887652 65211 33497 333755 015275 33871 30725 858698 822778 58511 41645 078888 014917 29864 72603 746788 837378 60887 15419 560673 814289 61093 18382 851329 851986 70019 58827 098884 813168 86384 49176 626898 863605 3680 26515 025453 811819 45319 61518 417312 074633 19064 68150 792614 810193 01198 17248 435837 083911 69718 22218 823395 888327 67415 76784 748725 091223 44282 51325 995868 086267 20483 34189 063570 096397 19272 77913 791268 004060 14298 00386 941331 099312 85991 85094 924786 001761 48871 39152 118312 n24 886485 68928 62605 626085 012793 81953 46752 156974 819111 88674 73616 389159 812583 74563 46828 296121 931504 26796 96163 374387 812167 04729 27883 120 043379 35076 26845 138487 0 11550 56688 53725 6 1353 054542 14713 88839 535658 910744 42701 15965 63 783 064809 36519 36975 569252 009761 86521 04113 888270 0 74912 41915 78554 364244 809619 01615 31953 275917 0 82088 19859 73902 921954 0 07334 64814 11080 3057 0 88641 55270 04401 034213 085929 85849 15436 780745 093827 45520 02732 758524 994427 74388 17419 806169 897472 85559 71309 498198 882853 13886 28933 663181 899518 72199 97021 360180 881234 12297 99987 199547 an 37833 85362 93550 43883 4224 67193 63625 13491 13443 70458 53363 78362 9036 01260 61574 1453 09996 15982 50581 08688 13403 86512 33 1EL41317311 5 PJlf dl lilCLAILITJIUEIXXIYPJ 70954771175351 075 75770144352 44704 Abscissasx7 Zeros of Laguerre Polynomials Weight Factorsw7 x 70 wiczi x w 0115031 quot7 2 719 1 336126 421798 039143 11243 16 1 411213 980424 092180 50285 29 1 199287 525371 148012 790994 2 474605 62765 208677 08055 3 559962 661079 277292 138971 3 6 6 8 4 1 3 hU NM 050570 54375 27 1053553 390593 153332 503312 01523 22277 341421 35523 73 1 145445 509407 445095 733505 o8g 0 r O h a U1 a J J 3994 Co COM WNW 713 1345523 59110 041577 45557 03 g1g711093 009929 107759 255927 1833359 899 0 N III 0 M 11076 933035 6 21227 541975 229420 03502 79 1 270517 733559 275214 295190 26 3740 18909 2 1 2908 403033 935321 023771 628994 50829 37 2103892 565016 560109 462543 3 032254 75096 19 1 503154 104342 083273 91230 30 7110 174575 11011 50 1 357410 592430 204010 243045 013779 34705 40 1 300441 115755 035400 97385 07 053662 0296 21 2 388879 085150 303114 630582 072945 45495 33 1K491119 929155 083190 23010 44 93950 09123 01 45394 705561 64814 508441 100034 29017 40 1 210050 207512 133020 055175 340143 35970 55 2 520074 550907 105300 390311 555249 51400 54 3 950151 597510 245025 555000 033015 27457 54 4 753000 300500 312275 415514 1104370 50379 00 5 202592 334950 393415 259555 5 wm000m700000 17 57909 40422 08 119 83956 482398 5 Eziiiig 23 g 8622 giiggi 363343 787350 2992059 70122 74 13 991102 721951 970459 504037 359542 57710 41 42 759424 495017 275944 324237 700501 00050 59 3 351175 057992 431555 590092 1254000 00442 75 5 233599 723050 721910 535435 nw12 44 90029020 15 99007 300010017 7 777 73 022204 55041 79 1 458964 873950 057353 55074 23 151251 02597 75 1 244002 011320 110770 139452 110893 21015 73 1 417000 830772 136925 259071 203375 13377 44 2 904492 222117 153045 423904 299273 53250 59 1 11333 382074 226068 459338 459922 70194 18 2 201023 811545 199032 750527 577514 35591 05 2 103991 974531 335052 458236 504452 54531 15 3 255397 354107 250074 575910 983745 74103 03 4 01017 202015 400502 500021 952131 50424 57 4 203231 592553 305532 151020 1590207 39305 02 89854 906430 784901 594560 1300505 49933 05 5 035505 505502 372320 911070 1711505 51074 02 7 1 55049 307554 4529 1 402990 2215109 03793 97 9 134239 103052 559725 045104 2840795 72509 84 12 305150 153504 721299 545093 7 3709912 10444 6 16 814807 74643 1054383 74619 n 019304 36755 60 1 400318 951701 04964 75975 40 102656 43953 39 1 421831 277862 117764 306086 255737 67449 51 1 147120 348658 191824 978166 495535 30845 25 2 205335 144507 277104 003523 amwww9wmm lmwwg mm 1273410 02917 90 5 8207 a 170 1 47 0 540543 243 009330 78120 17 1 218234 885940 02395 81703 11 193939533 78622 53 3 3 5 90 683 049259 17403 02 1 342210 177923 055010 00427 93 121559 54120 71 1 203027 577942 000700 02529 19 2 25994 95252 04 1 125425 810105 122355 440215 350752 27217 51 2 402050 549210 157444 072153 m8 i mi 017027 95323 05 1 369188 589342 043772 34104 93 1012022 05500 19 4 111574 392344 277404 192503 090370 17757 99 1 41886 780814 1 03386 93476 1313020 24021 70 5 545992 575202 325554 334540 225108 66298 66 1 15794 936537 166970 976566 1555440 77003 30 7 222531 590710 300531 171423 425570 01702 00 2 333434 922612 2 37692 47016 2077547 00994 49 9 422743 030490 445047 775304 704590 54023 93 3 79453 623523 320854 091335 2552309 42257 29 11 392109 725704 527001 770443 1075051 50101 01 5 907650 877336 42685 551083 3140751 91597 54 13 145551 520407 535955 345973 1574057 05412 70 7 8 4857 671627 5 81888 33635 3053050 3305 05 15 140302 705111 0 03173 7 3212 1 4 6 2205313 17350 09 9 104800 11987 690622 621529 4002500 55720 05 20 150059 490521 1152777 21009 Table 53 Abscissas and weight factors for Laguerre integration 51 1 SIN GULARITIES 070710 000000 122474 052464 165068 000000 095857 202018 043607 133584 235060 000000 081628 167355 265196 038118 115719 198165 293063 000000 xx 67811 00000 48713 76232 01238 00000 24646 28704 74119 90740 49736 00000 78828 16287 13568 69902 37124 67566 74202 00000 10187 32892 05845 32017 J39jwe 12fxdxiil w Abscissayix Zeros of Hermite Polynomials 86548 00000 91589 75290 85785 00000 13819 56086 27617 13697 74492 00000 58965 67471 35233 07322 46780 95843 57244 00000 52838 16668 31843 81528 107 n2 18 86226 92545 773 0118163 59006 1 295400 97515 n El804914 09000 2 813128 35447 775 19 45300 72040 4 93619 32515 2199532 42059 nz6 s lg724629 59522 1 157067 32032 3453000 99055 n 1 810264 61755 1 425607 25261 E Z 545155 82819 4971781 24509 1 661147 01255 1 2 07002 32581 2 170779 83007 4199604 07221 Table 54 04 09 55 25 29 N N 05 44 09 m o 68 HO wv t 95 82 bk HQ 14 61 wk J 26 392 wieti 146114 118163 132393 105996 124022 094530 098658 118148 087640 093558 113690 081026 082868 089718 110133 076454 079289 086675 107193 072023 073030 076460 084175 104700 Abscissas 11826 59006 11752 44828 58176 87204 09967 86255 13344 05576 83326 46175 73032 46002 07296 41286 00483 26065 01442 52156 24527 81250 27014 35809 034290 103661 175668 253273 343615 031424 094778 159768 227950 302063 388972 027348 024534 073747 123407 173853 225497 538748 1x 13272 08297 36492 16742 91188 03762 83912 26351 70805 70251 48978 10461 14491 85391 79909 21578 91619 79048 89393 07083 37285 62153 77121 40020 60584 3J3 om JBO ON Hm 24495 08900 and weight factors for If gxdxu ii 101 effgx Weight Factors 107 23705 89514 99882 32790 37738 54359 40164 52605 01060 20890 69782 O O 0 AW W Hermite 197 nr 10 1 610862 63373 1 2 40130 61108 2 338743 94455 3 134364 57467 67 64043 20552 570135 23626 260492 31026 516079 85615 390539 05846 857368 70435 265855 16843 7116 1 1 36455 integration 53 23 48 m 33 25 42 88 29 88 56 66 W DCD NHOU 0 U1 11 06 28 46 38 63 0 0 73 34 2 1007 068708 070329 074144 082066 102545 062930 063962 066266 070522 078664 098969 054737 055244 056321 058124 060973 065575 073824 093687 049092 049384 049992 050967 052408 054485 057526 062227 070433 089859 18539 63231 19319 61264 16913 78743 12320 27732 03661 39394 90470 52050 19573 78290 72754 69582 56728 56222 44928 15006 33852 08713 90271 03509 17423 24428 86961 29611 19614 35 513 049 Oh 5U 870 657 695 203 669 122 633 923 378 D05 BIN Mm 009 560 761 777 841 667 omoar u l HNampmla N buthOUIuH 769 532

### 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

#### "Knowing I can count on the Elite Notetaker in my class allows me to focus on what the professor is saying instead of just scribbling notes the whole time and falling behind."

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

#### "There's no way I would have passed my Organic Chemistry class this semester without the notes and study guides I got from StudySoup."

#### "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.