1Implicit, Implicit, MultistepMultistepand and Extrapolation MethodsExtrapolation MethodsLarry CarettoMechanical Engineering 501ABSeminar in Engineering AnalysisNovember 9, 20042Outline• Midterm solutions• Review last class– One step methods for numerical solution of differential equations– Local and global error• Multistep methods with constant and variable step size• Implicit methods using future time steps• Extrapolation methods3First Midterm Problem• Solve characteristic equation λ3+ λ2+ 4λ + 4 = 0 for λ = -1, ± 2iω0)0('',0)0(',1)0(0442233====+++ yyyydtdydtyddtydtCtBAeyt2cos2sin ++=−• Find A, B, and C from initial conditionstteyt2cos512sin5254++=−4First Midterm Problem, part b• Solution y = yH+ yPwhere yHis from part a with A, B, and C undetermined• Assume yP= D cos t + E sin t – Substitute this yPinto differential equation– Equate coefficients of sine and cosine on both sides of result– Solve for D and E0)0('',0)0(',1)0(cos442233====+++ yyytydtdydtyddtyd5First Midterm Problem, part b II• Solution y = yH+ yP= Ae-t+ B sin 2t + C cos 2t + D cos t + E sin t (know D and E)– Substitute solution into differential equation– Have three equations in three unknowns to solve for A, B and C– Result, not required for exam, shown below()tttteytcossin612cos1522sin154107++++=−6First Midterm Problem, Part c• Here the right-hand side is proportional to the sin 2t term in the homogenous solution• Use yP= D t sin 2t + E t cos 2t• Repeat process used for part b– Will get cancellation so that there will only be two kinds of functions in yPprocess0)0('',0)0(',1)0(2sin442233====+++ yyytydtdydtyddtyd27Second Midterm Problem• Solve using Laplace transforms0)0('',0)0(',1)0(2sin442233====+++ yyytydtdydtyddtyd[]42)(4)0()(4)0(')0()()0('')0(')0()(2223+=+−+−−+−−−ssYyssYysysYsysyyssYs()()()()()4144142)(2222++++++++=ssssssssY• Factor multiplying Y(s) same (with same roots) as λ equation from problem one4423+++sss8Second Midterm Problem II• Use partial fractions() ()415445111544125242524125104251011252)(22222222+++++++++−+++−+=ssssssssssssY• Use transform table and rearrange results for solutionttttttetyt2cos2532sin200932cos2012sin1012522)( ++−−=−9Review Numerical Approach• Solve initial value problem, dy/dx = f(x,y) (known) with y(x0) = y0– Use a finite difference grid: xi+1–xi= h– Replace derivative by finite-difference approximation: dy/dx ≈ (yi+1–yi) / (xi+1–xi) = (yi+1–yi) / h– Derive a formula to compute favgthe average value of f(x,y) between xiand xi+1– Replace dy/dx = f(x,y) by (yi+1–yi) / h = favg– Repeatedly compute yi+1= yi+ h favg10Review Notation and Order•xiis independent variable•yiis numerical solution at x = xi•fiis derivative found from xiyi: fi= f(xi, yi)•y(xi) is the exact value of y at x = xi•f(xi,y(xi)) is the exact derivative•e1= y(x1) – y1= local truncation error•Ej= y(xj) – yjis global truncation error• If e is O(hn), then E is O(hn-1)11Review Simple Methods• Huen’s method• Modified Euler method[]2),(),(),(2),(0111010111111101++++++++++++++=++=+=+=iiiiiiiiiiiiiiiiiiiiyxfhyyyxfyxfhyyhxxyxfhyy),(2),(2212111121121+++++++++=+=⎥⎦⎤⎢⎣⎡+=iiiiiiiiiiiiiyxfhyyhxxyxfhyy•Euler: yi+1= yi+ hifi= yi+ hi f(xi, yi)12Review 4thOrder Runge-Kutta• Uses four derivative evaluations per step),(2,22,2),(622311421131112111143211kyhxfhkkyhxfhkkyhxfhkyxfhkhxxkkkkyyiiiiiiiiiiiiiiiiiiii++=⎟⎠⎞⎜⎝⎛++=⎟⎠⎞⎜⎝⎛++==+=++++=++++++++++313Error versus Step Size for Simple ODE Solvers1.E- 161.E- 151.E- 141.E- 131.E- 121.E- 111.E- 101.E- 091.E- 081.E- 071.E- 061.E- 051.E- 041.E- 031.E- 021.E- 011.E+001.E-07 1.E- 06 1.E- 05 1. E-04 1.E-03 1.E-02 1.E-01Step size, hErroEul erHue nModifie d EulerRunge-Kutta0100=====−−xxatyywithedxdyyx14Error Propagation in Solutions of ODEs1.E-141.E-131.E-121.E-111.E-101.E-091.E-081.E-071.E-061.E-051.E-041.E-031.E-021.E-011.E+000 0.2 0.4 0.6 0.8 1xErroEuler (N = 10)Euler (N = 100)Euler N = 1,000)Huen (N =10)Huen (N = 100)Huen (N = 1,000)RK4 (H = 10)RK4 (N = 100)15Implicit Methods• Methods discussed previously are called explicit– Can find yn+1in terms of values at n– Use predictors to estimate y values between ynand yn+1• Implicit methods use fn+1in algorithm• Usually require approximate solution• Have better stability but require more work than explicit methods• Trapezoid method is an example16Derive Trapezoid Method• Subtract series expansion for yn+1about ynfrom series for ynabout yn+1)(2''321hOyhhfyynnnn+++=+)(2''31211hOyhhfyynnnn++−=+++())(2''''312111hOyyhhfhfyyyynnnnnnnn+−+++−=−++++17Derive Trapezoid Method II• Collect terms is last equation and substitute y’’n+1= y’’n+ hyn’’ + O(h2)())(2311hOhffyynnnn++=−++()())(4''''231211hOyyhhffyynnnnnn+−++=−+++()[])(4)('''''''232211hOhOhyyyhhffyynnnnnnn+−−−++=−++18Implicit Methods II• How can we use fn+1in algorithm to solve for the unknown yn+1?– One approach is trial-and-error solution– Euler step for first approximation of yn+1– Iterate on implicit methodnnnhfyy +=+)0(1()[]2)(1,1)1(1mnnnnmnyxffhyy++++++=419Implicit Methods III• Alternative approach for implicit method uses series expansion for fn+1about fn– Trapezoid method yn+1= yn+ h(fn+1+ fn)/2 + O(h3) example())(211hOyyyfhxfffnnnnnn+−∂∂+∂∂+=++()())()(2)(2321311hOhOyyyfhxfffhhOffhyynnnnnnnnnn+⎥⎦⎤⎢⎣⎡++∂∂+∂∂++=++=−+++20Implicit Methods IV• Neglect O(h3) terms and solve the last equation for yn+1-yn()⎥⎦⎤⎢⎣⎡−∂∂+∂∂++=−++ nnnnnnnnyyyfhxfffhyy112()22121hxfhfyfhyynnnnn∂∂+=⎟⎟⎠⎞⎜⎜⎝⎛∂∂−−+()nnnnnyfhhxfhfyy∂∂−∂∂+=−+21221• Have to compute partial derivatives of f(x,y) for implicit21Implicit Example• Look at sample equation dy/dx = f = -ay• Here, ∂f/∂x = 0 and ∂f/∂y = -a() ()hayhahahayhayahhayyyfhhxfhfyynnnnnnnnnn+−=+−+=−−+−+=∂∂−∂∂++=+22222)(2022221• Here yn+1= G ynwith G = (2 – ha)/(2 + ha)• Use next week to discuss stability22Another Implicit Approach• Use Newton-Raphson iteration for yn+1– Solve g(x) = 0 by iteration x(m+1)= x(m)–g(x(m)) / g’(x(m)) –g(yn+1) = yn+1–yn–hfn/2 – hf(xn+1,yn+1)/2– g’(yn+1) = fn+1–0 –0 –h(∂f/∂y)/2 ()())(1)(11)(11)(1)(1)1(12,2,2mnmnnmnnnnmnmnmnyfhyxfyxhfhfyyyy+++++++++⎟⎟⎠⎞⎜⎜⎝⎛∂∂−−−−−=23Multistep Methods• Previous methods used only information