COMP METH CIVIL ENGR
COMP METH CIVIL ENGR CGN 3421
Popular in Course
Popular in Civil Engineering
This 34 page Class Notes was uploaded by Geraldine Kuhic on Friday September 18, 2015. The Class Notes belongs to CGN 3421 at University of Florida taught by Staff in Fall. Since its upload, it has received 17 views. For similar materials see /class/206921/cgn-3421-university-of-florida in Civil Engineering at University of Florida.
Reviews for COMP METH CIVIL ENGR
Report this Material
What is Karma?
Karma is the currency of StudySoup.
Date Created: 09/18/15
mum Emmxnh mh any Numerical Methods Lecture 6 7 Optimization Tapas numznulapnmuman Newmnagam Randnm gm GaldaneconnSearch Op mizzvim r mummuun Wm unsungwhznsamzfmmnnuhzsamaxxmnmmmxmmum mdxwhzxe k max ax x 7 mm 3 Faxexam z A mm xepxesemsth castafapmject mmnmu m cast afthz mm thm thnanexactsalmmms mummy nub pun 2mm 1 Gwen y 56X73j2Fmdxvhznyxsmxmmnm mlyncalmhmnn 1 12x733 u Ekamgltzl Gwen y Emmnwizsmo swinmmnymmm Defends anthz range we39re mmam Hm 1mm maximandmmmamzms havmg 1mg afzzmslnp cases Anexactsalmmn wuldbe a s m pain 7 a u s 39u mmmmms mum m Imam CGN 3421 Computer Methods Gurley Single variable Newton Recall the Newton method for nding a root of an equation Wk d x xk1xk where f xk 2 16 f xk We can use a similar approach to nd a min or max offx The min max occurs where the slope is zero So if we nd the root of the derivative we nd the max min location f 00 gm o Find roots of gx using Newton 80 xkw x x k1 k v x g k f 06 Example Newton Method mumquot 3C 1 y V g to find where the nd the maximum of this function 239 X occurs N 4 3 2 fx 6 5x 2x 2416 1U we ll need these B 7 7 D D5 3911 15 2 25 393 f x 4x 15x 4x24 w 9 1 it 2 EDT find the root of the f X 1216 30264 function s derivative for x 03 G 25 I I I I OR B 05 u 1 5 2 25 3 we can use central difference to replace the analytical derivatives Numerical Methods Lecture 6 Optimization page 104 of 1 ll CGN 3421 Computer Methods Gurley single variable Random search A brute force method 1 Sample the function at many random x values in the range of interest 2U Q n a I n p H 5 2 If a su icient number of samples are selected 15 iteration 2 a number close to the max and min will be 39 found 10 quot 3 Reset the range to a smaller subrange and look again Now the values are more dense and you are more likely to nd the max 35 1 L5 2 A simple Mathcad program for random search ORIGIN E 1 num 10 fx x4 5x3 Zx2 24x x runifnum 03 fx fx biggest maxfx FindiMaxilndex fx index lt l for i E 2 length fX index FindiMaxilndex fx index lt i if fx gt fx 1 index index Start over and narrow region x runifnum xmdeX 25xmdeX 25 fx fx biggest maxfx 2 a 1L 1 639 index FindiMaxilndex fx i l Note that if we make num bigger the rst and second iterations will be much closer since the density of random numbers over the same area will increase improving the liklihood of hitting the largest number Numerical Methods Lecture 6 Optimization page 105 of 1 ll mum Emmxnh mh mm mm nanafthz cad xsl mm Note that we39ve zoomed m dvanrgtz ozyandom march 51mph m xm zmzm 015mgme glnbal rmm lacslmaxxmals easy mm m equtmnwvexalngexmnge Fm mummm Newmn leaves wnhlmsafmmma andmaxxma1fNe ncmmnf mm Randnm Search can simplyplck dqmth and D m mam maximum Smce Itcampnresalhhz xwlzs anyway 5 o 5 quot3 u 75 o 5 m Dtmdwznmgtz ogmdom Search am farce methndr Ises m span mfannmmn sham m mm Dr 5 dzths mnnmhexpummmymlssmnaw m Becamzsvexyslnwfaxm nrmmznsmmlpmblzms m get a gaad 1m 1 fp q x y med mun4 nnmhexsm gen gaadlnak fm medsm nnmhexs mmmmms mum m m vi 111 CGN 3421 Computer Methods Gurley Single Variable Golden Section Search Optimization Method Similar to the bisection method Define an interval with a single answer unique maximum inside the range sign of the curvature does not change in the given range Divide interval into 3 sections by adding two internal points between ends Evaluate the function at the two internal points X1 and X2 if fX1 gt fX2 the maximum is between an and X2 rede ne range an an XmX X2 if fX1 lt fX2 the maXimum is between X1 and XmX rede ne range an X1 XmX XmX I I xmn X1 X2 me xmn X1 X2 me look inside red range for second iteration Divide new smaller interval into 3 sections and start over Q1 Where should we place the two internal points Let s see ifwe can pick some efficient points A L xmn X1 X2 me 4 L1 L2 Set the following conditions 1 L L1 L2 2 R LL2 L2L1 substitute 1 into 2 gt 1 R lR solve for R R is called the Golden Ratio named by ancient Greeks R sqrt512 61803 So ifthe range is 0 1 X11 R 38197 Numerical Methods Lecture 6 Optimization page 107 of 111 CGN 3421 Computer Methods Gurley X2 0 R 61803 3 8 197 6 l 803 0 l i i i i For range 0 1 The ratio of the larger length to the smaller remains L2 Ll R 61803 So in general for the range is xmn xmx Xl xmxR xmxxmn X2 xmn R xmx xmn xmn xmx Rxmx xmn xmn R xmx xmn xmx i i For generlc range xmn xmx Why is this particular choice of internal points ef cient Because we enforce the ratio LL2 L2Ll every successive iteration reuses one of the previous internal values Q2 How do we evaluate when we are close enough to maximum to stop We can t use err lfxl like we did with root nding algorithms since the value of fx says nothing about whether its a maximum value We ll instead evaluate error with respect to the xaxis xopt xopt m5 ERR 39 OR gt lt gt range range I xmn x1 x2 xmx xmn x1 x2 xmx The maximum x distance between the current guess and the current range is called the error ERR lR xmx xmn Numerical Methods Lecture 6 Optimization page 108 of 111 CGN 3421 Computer Methods Gurley So ifwe set the tolerance to 001 the nal X guess is within 001 of the X location of the actual maximum A few iterations Km X1 X2 XmX an X1 X2 XmX an X1 X2 XmX X1X2 an XmX Numerical Methods Lecture 6 Optimization page 109 of 1 ll CGN 3421 Computer Methods Gurley A working algorithm GoldenSearch f a b tol grat 52 1 d grat b a X2 a 1 X1 b d err 100 numit 0 while err gt tol numit numit 1 if f X1 gt f X2 Xopt X1 err 1 grat b a Xopt if err gtt01 b X2 X2 X1 1 grat b a X1 b 1 otherwise Xopt X2 b 3 Xopt err 1 grat if err gtt01 alt X2 X1lt X2 dlt gratb a X2lt ad Xopt j f Xopt Numerical Methods Lecture 6 Optimization page 110 of 111 CGN 3421 Computer Methods Gurley 2 fX 2sinX out 2 G01denSearchf04 01 El Reference C 11ine11athcadTut0ria1 sMyFuncti0ns mcd XI CreateVect0r0 4 1 fX outz OO X outl Numerical Methods Lecture 6 Optimization page 111 of 111 CGN 3421 Computer Methods Gurley Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods Lecture covers two things 1 Solving systems of linear equations symbolically 2 Using Mathcad to solve systems of nonlinear equations 3 Investigating algorithms to find roots of equations Solving Systems of Linear Equations Symbolically Let s take a look at a very powerful tool in Mathcad that started a revolution in computational analysis It started in the late 1980 s when I was an undergraduate A company called Wolfram created a computer program called Mathematica This was the rst computer code that could solve algebraic and calculus equations symbolically That is if I had an equation that said xy z Mathematica could tell me that y z x without ever needing me to assign numbers to x y or 2 It also was able to solve integrals dif ferential equations and derivatives symbolically This was an incredible advance and opened the doors to a whole new world of programming numerical methods pure mathematics engineering and science Since then a competing code called Maple was developed and sold itself to other software companies to include in their programs The end result Mathcad uses Maple as a solving engine in the background you don t see it to solve problems symbolically Here we will look at a brief example of how to use this capability in the context of solving a system of linear equations Example The structural system below is something you will see in CBS 3102 or CES 4141 r1 and r2 are labels that indicate how the ends of the beam are allowed to move Q and W represent the external loads a couple and a distributed load respectively and material properties are given as E and I L is the length of the beam The goal is to solve for the amount or rotation at rl and the de ection at r2 that occurs for thegiven loads This would help us to solve for internal stresses allowing us to design the beam to sur vive these internal forces W3KFT w Q Q2KFT C L15Fr E4000ksi L L I 1400 inA4 r1 1 r2 Solution The way we learn to solve this problem in CES 4141 is using a Matrixbased solution proce dure The generic form of the solution is K r R K is a 2 x 2 stiffness matrix that contains information about the structures shape boundary conditions and material properties THe information needed is L E and I r is a 2 x 1 vector that contains the unknown rotation and displacement quantities sought R is a 2 x 1 vector that contains only information about the external loads W and Q for this problem Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 68 of 82 CGN 3421 Computer Methods Gurley Solution continued So we will m K as a known 2 X 2 stiffness matrix We will m R as a known 2 X 1 load vector We will m for the unknown displacement vector r We will not go into any detail on HOW we ll in the values for K and R that s for another class Here we will just consider how to solve the ssytem of matrix equations symbolically Define Load Vector R as a Define Stiffness K as a function of E L function of Q W L 4131 6131 known information wL2 L L2 Q T KEIL R WL 6EI 12EI Q wL L2 L3 20 The variables in the parenthesis EL and QWL is telling mathcad that the variables K and R respectively will be a function of those variables Show me the inverse of the stiffness matrix Not a necessary step but pretty cool see the matrix version of 1K Symbolic Inverse 1 1 2 Note the arrow gt is the way to E1 39 2E1 39L initiate a symbolic operation Here KE1L 1 a we are asking what is the inverse of 1 L2 1 L3 K The arrow comes from the 2131 3131 view gt toolbars gt symbolic pad Calculate the displacement vector as K391 R Note that we DO need to keep writing the names of the independent variables next to K r and R to tell Mathcad what symbols to look for 1 S bolic Solution rEIL W KEIL R wL 3quot Q Q Symbolic Calculation Display the resulting displacement vector calculated in terms of W Q L E Note the arrow displays a symbolic result Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 69 of 82 CGN 3421 Computer Methods Gurley We can now send specific values into the function for r by first assigning numbers to the parameters Note that units are not being declared so its up to the user to be consistent Here are the mate rialgeometric E 4000 1 1400 L 1512 properties W i Q 212 Here are the external load values 12 Finally we can now ask Mathcad to evaluate the results of sending the parameters into the function we created above Now that E l L Q W have numbers they will be subsituted into the equations to get answers Be sure to use the same order for the parameters as defined above w Q412 Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 70 of 82 CGN 3421 Computer Methods Gurley Solving Systems of Nonlinear Equations We won t go into the algorithms themselves here We will just focus on how to use Mathcad to solve the problem Use Mathcad help and use the keywords nonlinear equations to get some information Follow the link to the Find function in Mathcad Find is the workhorse that gets the solution All you need to do is give it the correct description of the problem This includes the equations to solve any possible constraints to the solutions you are interested in and a place to start looking for a solution ie initial guess for the solution Here we will go straight into an example to show how it s done ORIGIN 1 Using the Find function to solve systems of nonlinear equations This example will solve for the intersecting values of the followin system of 2 equations Let39s first take a look at the solution visually by plotting We will create some functions that allow us to plot the two equations above I 2 2 xy2 839 f1xquot87xz2 X 15X2y1 39 2 574994 f2 8 X 2 i f3x 1 7x7 15x2 i A u 1 10 E Fa can see that there are two solutions to the two equations above We39ll try to find both solutions by using the 39find39 function Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 71 of 82 CGN 3421 Computer Methods Gurley Use the help desk with the keyword 39find function39 Step one start by guessing the solution to the set of two equations we are trying to solve Let39s guess that x4 and y 1 Next we create what is called a solve block This defines the equations to be solved 1 type the word 39Given39 2 type below 39Given39 all the equations to be solved Note that we use the 39boolean39 tab to get the equal sign not the keyboard stroke Given X2 y 7 22 8 X 1539X2 y 1 Finally we use the 39Find39 function to solve by sending in the two guesses we made above for x and y 32 1 i v mm Am 39Find39 starts looking at the initial guess values and iteratively updates the values ofx and y until a pair is found that solves the equations defined in the solve block The results give the x and y values corresponding to one of the solutions find where x1278 and y 0523 in the graph above 39Find found this solution because the initial guesses were closer to this solution Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 72 of 82 CGN 3421 Computer Methods Gurley Now let39s find the other solution seen in the graph we39ll redifine the initial guesses to somewhere near the other solution and use 39Find39 again X2y 228 X15X2y1 Notice now that the solution changes to the other one shown in the earlier graph What else can we do with this solve block Suppose we want only to nd solutions ifx gt O that is we want to tell Mathcad to not even bother with solutions ifx lt 0 We can do this by adding constraintsto the solve block X 4 yz2 Given X2y 228 X15X2y1 Xgt 0 Adding the constraint X gt 0 i373 Solution changes to the first one we found earlier even though the initial guess is closer to the solution with the negative X value a vquot L LEV rm 5 J y e Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 73 of 82 CGN 3421 Computer Methods Gurley Solving for Roots of Equations What Locating where some function crosses the XaXis Find X such that 0 Why We can then find where crosses any value Find X such that 4935 or Find X such that gx 4935 0 When Numerical methods in general are 1 less accurate 2 slower than using an analytical eXact solution Use numerical methods when analytical answers not available EXample Given y 5 626 Find X such that y 0 analytical solution 0 5 626 gt x 56 no need for Given y exp sinx 3 sin15x 5 Find X such that y 0 No idea Good a prime candidate for a numerical solution We will consider several methods each does the following in some way y error y distance from 0 1 Start with some initial guess or range in I searchinthe range lb for afoot which to look I b 1 rst guess 2 Use some information about the function x to make another closer guess 3 Stop when our current guess is close enough to a solution new error y distance from o I b a X The user will decide what is close enough next guess how big the error can be before we can stop looking Method 1 INCREMENTAL SEARCH METHOD Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 74 of 82 CGN 3421 Computer Methods Gurley 1 Start with an initial range that contains the root and subdivided it into several smaller subranges 2 Look inside each subrange one by one for the root When the subrange containing the root is identi fied choose either end ofthe range as the guess 3 Evaluate the error in your guess If error is too big subdivide your new range into smaller subranges 4 Loop back to step 1 and continue the loop until the error in step 3 is small enough start looking in range a1 b1 segment into 4 pieces we picked 4 out of thin air subrange b2 a1 a2 subra subrange b1 subrange subrange b3 b2 a2 a3 subra subrange subrange look in each subrange until a2 b2 is found evaluate the error err minimum l fa2 l l fb2 l Now divide a2 b2 into 4 segments and look again look in each subrange until a3 b3 is found evaluate the error err minimum l fa3 l l fb3 l Now divide a3 b3 into 4 segments and look again Out of the 4 subranges we pick the one in which the function crosses the XaXis as the new range to subdivide even smaller Visually we can see this but how does the algorithm know how to identify the smaller range Remember we have the function so we can evaluate it at any X location we choose Take a moment work out a way to identify the new smaller subrange Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 75 of 82 CGN 3421 Computer Methods Gurley 8 2 subrange 3 3 5 4 2 2 5 subrange subrange subrange 6 8 2 subrange 3 3 5 4 Z 2 5 subrange subrange 6 Incremental Search Root Finding Algorithm ncSearchnumsegs tolfuncab err lt mirfunca Ifuncb funcX 2 5 2X2 while er thol i2 1 I9N3 dx 9 m numsegs for i E l numsegs xlt a if ldx prod lt funcxfuncx dx if prod S 0 10 funci 0 ex b lt X dx err e mirfunca Ifuncb i root lt a if Ifunca lt Ifuncb root lt b otherwise root Something new to notice in how we use the above functionwe ll discuss in class the idea of passing a function into another function Note that the equation we are solving in the above example is not explic itly typed intop the IncSearch function Is this a good thing Why This is a brute force method that uses very little information about the function to improve the guess Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 76 of 82 CGN 3421 Computer Methods Method 2 BISECTION METHOD 1 Start with an initial range a b and cut the range in half assume this half point m is the root 2 Evaluate the error err l fm l 3 If error is too big decide if the root is to the left or right of m If root is to the left ofm reassign a new rangeaa bm Else Ifroot is to the right of In reassign a new range a In b b 4 Go back to step 1 stop bisecting when the error at m is small enough Three Bisection iterations Eventually we zero in on the root How do we decide where the root is relative to m for the next a b Same as the incremental search prod fa fm if prod lt 0 root is to the left otherwise root is to the right Gurley m b err lfml err lfml new b new In evaluate the error err l fm l err gt tol yes make new range root is to left ofm gt a a b m root is to right ofm gt a In b b divide new a b in half to get new m evaluate the error err l fm l etc iteration 1 7 m M error at ml too big so guess again fml iteration 2 b2 error at m too big so guess again iteration 3 N 13 b3 33 Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 77 of 82 CGN 3421 Computer Methods Gurley This slow method does not use much function information We can do better Method 3 FALSE POSITION Relative height of function at end points used to make better guesses fa Given the heights of fa vs fb it seems more likely that the root would be closer to a than b and not in the middle Let s use this concept in an algorithm fb 1 Define initial range a b possibly the result of a single pass of the incremental search method rst guess c intersection of a line from fa to fb with XaXis c a b b a from similar triangles f b fa 2 Evaluate the error err l fc l If error is too big decide if the root is to the left or right of c for fa fc lt 0 reassign a new range a a b c for fa fc gt 0 reassign a new range a c b b new guess 6 a b b a f b fa 3 Stop iterating when the error is small enough Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 78 of 82 CGN 3421 Computer Methods Gurley Two False Position iterations iterationl error at m2 too big Eventually we zero 1n on the root c180 guess again b1 Tends to zero in faster than using bisection method why 31 This method uses more function information el than biSeCtion a stays b moves over from right bisection only knows if fa fm are positive or negative 02 iterati on 2 false position Uses of fa and fc and uses their relative magitudes to make next 32 guess at c Method 3 Root nding Newton Raphson method Incremental search uses sign of fa and fb Bisection uses sign of fa and fb False position uses sign and relative magnitude of fa vs fb Newton Raphson uses derivative of fx dfxdx error y distance from 0 1 Pick a single point not a range as an initial guess xl tangent to xl second guess X 2 evaluate the error err fxl 2 X1 rst guess l 3 Next guess draw tangent line from initial guess to the xaxis New guess x2 is the intersection of the tangent line with the xaxis 4 Back to step 2 to evaluate error of new guess etc Tangent line slope is given by derivative of fx General formula for next guess xkl in terms of previous guess xk fltxk newguess Xk1 Xk W Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 79 of 82 CGN 3421 Computer Methods Gurley this is from the de nition of slope fXk tangent to fXk 0 fxk slope slope slope dfxk dquot quotk 1 quot k solve for Xk1 Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 80 of 82 CGN 3421 Computer Methods Gurley Let s continue with a few more iterations error y distance from 0 tangent to fX1 second guess X2 X1 rst guess l second guess I X I 2 X1 X3 tangent to fX2 I X3 1 X4 X2 3 X1 X Another example 2 iterations in one gure slope f X2 slope f Xl X1 rod Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 81 of 82 CGN 3421 Computer Methods Gurley Some Comments 1 Quickly converges to the root under the right conditions 2 Can be divergent a very bad word if initial guess not close to the root Must have condition in the inde nite loop to stop if divergent may need incremental search to narrow the range example X1 is not close to the root tangent line send us off in wrong direction we re not getting closer f X2 3 Requires the function AND the function s derivative Subtle Foreshadowing What if an equation for the derivative dfX is not readily available Sure would be nice to have an algorithm that numerically estimates the derivative of any function HHmmmm Good test question Under what conditions will Newton Raphson nd the exact root after a single iteration Numerical Methods Lecture 3 Nonlinear Equations and Root Finding Methods page 82 of 82 CGN 3421 Computer Methods Gurley Numerlcal Methods Lecture 4 Numerlcal Integratlon leferentlatlon Numerical Integration approximating the area under a function We will consider 3 methods Trapezoidal rule Simpson s rule Gauss Quadrature The basic concept is to break the area up into smaller pieces with simple shapes tted to approximate the func tion over short lengths The area under these simple shapes are then added up to approximate the total Trapezoidal rule piecewise linear approximation to curve fx A112dxfx1fx2 chop area into N smaller pieces sample of two pieces on the right connect areas with straight lines creating trapezoids nd the area under each trapezoid A12dxfxfxdx add up the trapezoids A1 Example the gure at the bottom of the page chops a function into 4 trapezoids a X f X over the range a b The trapeziodal rule adds up the area of each of the 4 trape zoids to represent the area under the fucntion fx x1 X2 X3 dx baN a b are limits N is number of pieces we divide the area into A1 12 dx fx1 fx2 A2 12 dx fx2 fx3 A3 12 dx fx3 fx4 A4 12 dx fx4 fx5 Area 12 dx fx12fx22fx32fx4fx5 N The generalized equation expressing trapezoidal rule is then gt Area 5 dxfx1 2 Z xi fxN1 i 2 Error The error in blue between the area of the trapezoid and the curve it represents builds up Error reduction We can reduce error if we can shorten the straight lines used to estimate curves Increasing the number increase N 0ftrapez0ids used over the same range ab will reduce error in the area estimate 5 Numerically estimate area under fx within a b Numerical Methods Lecture 4 Numerical Integration Differentiation page 83 of 88 CGN 3421 Computer Methods Gurley We can also get a little more sophisticated than straight line estimates The error using trapezoidal rule comes from using a series of straight lines to represent a curve To reduce this error we ll now look at using a series of second order curves to represt the function Simpson sRule piecewise second order approximation of fX AZdX3fXi14fXifXi1 chop area 1nto N smaller p1eces connect every 3 points with second order curve fit this gives us N2 smaller areas representing the curve in the range ab Find the area under each second order polynomial estimate Adx 3 fX14fX2fX3 Add up the smaller areas to estimate area under the curve Xll X1 Xl l l Example the gure at the bottom of the page chops a function into 8 pieces over the range a b Every two adjacent pieces are joined to form a subarea The Simpson s rule adds up the 4 subareas to represent the area under the fucntion fX A1 dx3 fX14fX2fX3 A2 dx3 fX34fX4fX5 A3 dx3 fx54fx6fx7 A4 dx3 fx74fx8fx9 Area A1A2A3A4 dX3 fx14fx22fx34fx42fx54fx62fx74f8fx9 N N l 1 Slmpson s rule 1s gtArea E gdxfx1 4 Z fxJ 2 Z fxk fxN1 j 2 even k 3 odd Reduce error by increasing N This method requires N to be even since we need 3 points per subarea fX Numerically estimate area under fX within a b N 8 segments gives us 4 areas to add up each area is a topped off with a second order polynomial red line a l 1X2X X4 Numerical Methods Lecture 4 Numerical Integration Differentiation page 84 of 88 CGN 3421 Computer Methods Gurley Gauss Quadrature Integration Rather than cutting the range up into small pieces we can try to nd a single trapezoid that represents the function over the given range The trapezoid is defined by the range ll and the evaluation of the function at fxl and fx2 See figure gt If we pick x1 and x2 correctly the blue area under estimate will balance out the red area overestimate and the total area under the trapezoid will be the same as the function r 1 1 The trick is then to nd the correct values for x1 and x2 and to weight the function evaluations Our integral estimate becomes 1 I 1fxdx cofx0clfx1 We need to solve for 4 unknown values CO Cl x0 x1 We ll need four equations We ll get these by assuming that the above estimate is exact for constant through 3rd order functions 1 l cofx0clfx1 11dx 2 cofx0clfx1 J1xdx 0 1 2 2 1 3 cofx0clfx1 I 1x dx 3 cofx0clfx1 1 1x dx 0 solving these four equations simultaneously gives us the following solutions 1 1 1 CO Cl x0 261 J3 The result This says that we can estimate the area under any function between l l Just by evaluating the fiunction at two carefully chosen spots 1 1 1 mm f f lt1 I 1 J3 J3 Also the derivation provides that the answer is exact if the function being integrated is a polynomial up to third order The extension What if the range of integration is not l l Let s say 6 2 when the range is l 1 Let s also say the range of interest is a b Skipping the derivation We reduce the integration to an equivalent over the range l l by applying the following change of variables Numerical Methods Lecture 4 Numerical Integration Differentiation page 85 of 88 CGN 3421 Computer Methods Gurley x W 2 dx bid 3 2 To apply this change of variables we rewrite equation 1 now as a i a a i xw f 29 m f 29 m 3 ZZ4 m Integrate the following function over the range 0 08 8 IO 02 25x 200x2 675x3 900x4 400x5 dx 0 Rewriting Eq 4 above gives the general solution ftde fx0fx1b a 5 For this example a 0 and b 08 Plugging these values into 6 then plugging x0 x1 into 5 we get Gauss Quadrature Area estimate f0l69 f063l 0516741l305837 182257 Exact Answer 1640533 3 Note if we truncated the function be1ng 1ntegrated after the 6 term then the area we estimate usmg 5 and 6 will be exact since Gauss quadrature is exact for up to third order polynomial equations Extension Gauss quadrature can be extended to higher orders involving more evaluation points for example 1 I 1fxdx cofx0 clfx1 cz xz cn 1fxn 1 Numerical Methods Lecture 4 Numerical Integration Differentiation page 86 of 88 CGN 3421 Computer Methods Gurley Numerical Differentiation Read section 231 Derivative represents the slope of a function fX at a particular value x slope is Lse change in fXchange in X run formally the de nition of the derivative is f x lim fXAX ZX lim fXAX ZX AX Ax gt 0 A36 Ax gt 0 ZAX Here we will approximate the limit by taking some small Ax rather than mathematically letting it go to 0 Forward dz erence method f X liIII f MAX x Ax gt 0 A26 for ward difference nd the derivativ at X X Backward difference method XdX M lim fltxgt ltx Axgt b Ax gt 0 A36 nd the derivative at X ackward dif erence XdX X Central di erence method lim fXAx fX AX Ax gt 0 ZAx Central dlfference nd the derivative at X XdX X XdX Numerical Methods Lecture 4 Numerical Integration Differentiation page 87 of 88 CGN 3421 Computer Methods Gurley Comparison of all three methods backward differenc central difference forward difference nd the derivative at X XIdX E xdX Typically central difference is the most accurate method since forward and backward biases are averaged out In fact just like Simpson s rule is one order higher accuracy that Trapezoidal rule the central differ ence method has one order higher accuracy that forward or backward difference User supplies location X at which to estimate derivative User supplies the Chi value AX Apply one of the three equations to estimate derivative Mixing and Matching Numerical Methods Consider the Newton Raphson root finding method from NM Lecture 3 This requires the function and its derivative to iteratively solve for a root Suppose Ihave a complicated looking function and I don t have the time skill or inclination to find its derivative analytically to use in the Newton Raphson routine How might I combine the central difference method with Newton Raphson Numerical Methods Lecture 4 Numerical Integration Differentiation page 88 of 88 CGN 3421 Computer Methods Gurley Lecture 4 Example Algorithm Development Another 2 algorithms swapping and sorting We saw at the end of Lecture 3 how to develop an algorithm to identify the maximum value out of a vec tor and the location of that value In this lecture we will combine that idea with the ability to swap values within an array to create a simple sorting algorithm Remember that there are already builtin Mathcad commands to sort a vector The purpose of this exercise is to learn the basic use of control structures and algorithm development so that we can create our own more advanced programs later on swapping switch the contents of two variables 3 Given x 3 but I want x 2 ie I want to swap the contents of spots 1 and 2 1 1 How can I swap the contents of X1 with X2 Will this work x1 x2 x2 x1 Why not x12 x2 gt xL 3 x2 xL gt x2 3 What we need is some place to temporarily store one or the other temp x1 x1 x2 x2 temp 2 temp x1 gt temp2 x 3 1 3 x1x2 gt temp2 x3 1 3 x2 temp gt temp2 x 2 1 mission accomplished what would pseudo code look like Lecture 4 Example Algorithm Development page 39 of 46 CGN 3421 Computer Methods Gurley 1 identify the variables to swap 2 save the contents of l to a temporary location 3 replace contents of l with that of 2 4 replace the contents of 2 with that in the temporary location From now on this entire pseudocode can be a one word command SWAP variables 1 and 2 selection sort reorder contents of a vector from low to high or high to low Combines the previous two algorithms 1 nding a maximum 2 swapping 41 we have student grades G 17 we want to reorder from LOWEST to HIGHEST Pseudo code 1 Find the location of the lowest value in the entire vector 2 Swap the contents of that spot with the 1 spot 3 nd the location of the second lowest value in the entire vector a or nd the lowest value in the entire vector excluding 1 4 Swap the contents of that spot with the 2 spot 5 Find the location of the lowest value remaining in the entire vector excluding spots 1 and 2 6 Swap the contents of that spot with 3 7 repeat steps 5 and 6 moving down one spot each time Sorting let s illustrate how we intend this algorithm to work 41 73 Given G 17 52 13 13 13 1 3 73 17 17 after iteration l G 17 after iteration 2 G 73 after iteration 3 G 41 5 2 52 5 2 4 1 41 7 3 Lecture 4 Example Algorithm Development page 40 of 46 CGN 3421 Computer Methods Gurley after iteration 4 G 41 we re done of iterations needed One less than the number of values in the vector This information will be used when we set up a loop Making it Work Let s build this up slowly by re ning the pseudocode and putting together one iteration at a time Pseudocode for iteration 1 only Fill in some Mathcad code for iteration 1 Repeat for iteration 2 Repeat for iteration 3 Generalize when we find a pattern Lecture 4 Example Algorithm Development page 41 of 46