1Introduction to Numerical Introduction to Numerical SolutionSolutionof Ordinary Differential Equationsof Ordinary Differential EquationsLarry CarettoMechanical Engineering 501ABSeminar in Engineering AnalysisOctober 28, 20042Outline• Review last class and homework• Midterm Exam November 4 covers material on differential equations and Laplace transforms (no phase plots)• Overview of numerical solutions– Initial value problems in first-order equations– Systems of first order equations and initial value problems in higher order equations– Boundary value problems– Stiff systems and eigenvalues3Review Last Class• Phase plots, critical points, and stability• Look at system of two linear homogenous, autonomous equations–dy/dt = Ay (no function of time)• Critical points and stability depend on matrix eigenvalues which depend on determinant properties• Describe node, center, saddle point and spiral4Numerical Analysis Problems• Numerical solution of algebraic equations and eigenvalue problems• Solution of one or more nonlinear algebraic equations f(x) = 0• Linear and nonlinear optimization• Constructing interpolating polynomials• Numerical quadrature• Numerical differentiation• Numerical differential equations5Interpolation• Start with N data pairs xi, yi• Find a function (polynomial) that can be used for interpolation• Basic rule: the interpolation polynomial must fit all points exactly• Denote the polynomial as p(x)• The basic rule is that p(xi) = yi• Many different forms6Newton Polynomials• p(x) = a0+ a1(x – x0) + a2(x – x0)(x – x1) + a3(x – x0)(x – x1)(x – x2) + …+ an-1(x – x0)(x – x1)(x – x2) … (x – xn-2)• Terms with factors of x – xiare zero when x = xi– Use this and rule that p(xi) = yito find ai•a0= y0, a1= (y1–y0) / (x1–x0) •y2= a0+ a1(x2–x0) + a2(x2–x0)(x2–x1) – Solve for a2using results for a0and a127Newton Polynomials II•y2= a0+ a1(x2–x0) + a2(x2–x0)(x2–x1) ))(()())(()(1202020101021202021022xxxxxxxxyyyyxxxxxxaaya−−−−−−−=−−−−−=• Could continue in this fashion to determine coefficients from data • Use alternative scheme – not derived here – known as divided difference table to compute akfrom same data8Divided Difference Table• Enter data on xiand yiin rows of table skipping one row between entries• Start with yidata as zeroth divided difference• First divided difference, Fi= (yi+1–yi) / (xi+1–xi)– Second (or later) divided difference is difference of first (or later) differences–aicoefficients are initial divided differences9Divided Difference Tabley3x3↑ a3y2x2←a2y1x1←a1←a0y0x001010xxyyF−−=12121xxyyF−−=23232xxyyF−−=02010xxFFS−−=13121xxFFS−−=03010xxSST−−=10Divided Difference Example10030↑ a34020←a21010←a1←a00010100100=−−=F3102010401=−−=F62030401002=−−=F1.020130=−−=S15.1030361=−−=S60010301.15.0=−−=T11Divided Difference Example II• Divided difference table gives a0= 0, a1= 1, a2= .1, and a3= 1/600• Polynomial p(x) = a0+ a1(x – x0) + a2(x –x0)(x – x1) + a3(x – x0)(x – x1)(x – x2) = 0 + 1(x – 0) + 0.1(x – 0)(x – 10) + (1/600)(x – 0)(x – 10)(x – 20) = x + 0.1x(x – 10) + (1/600)x(x – 10)(x – 20)• Check p(30) = 30 + .1(30)(20) + (1/600) (30)(20)(10) = 30 + 60 + 10 = 100 (correct)12Divided Difference Codefor ( i = 0; i < n; i++ )D[0][i] = y[i];for ( k = 1; k < n; k++ )for ( i = 0; i < n – k; i++ )D[k][i] = ( D[k-1][i+1] –D[k-1][i] ) / ( x[i+k] – x[i] );• D[k][i] is ithvalue of kthdivided difference• Code for n data points (0 to n–1)313Constant Step Size• Divided differences work for equal or unequal step size in x• If ∆x = h is a constant we have simpler results–Fk= ∆yk/h = (yk+1–yk)/h–Sk= ∆2yk/2h2=(yk+2–2yk-1+ yk)/2h2–Tk= ∆3yk/6h3= (yk+3–3yk+2 + 3yk+1 –yk)/6h3– ∆nykis called the nthforward difference– Can also define backwards and central differences14Interpolation Approaches• When we have N data points how do we interpolate among them?– Order N-1 polynomial not good choice– Use piecewise polynomials of lower order (linear or quadratic)– Can match first and or higher derivatives where piecewise polynomials join– Cubic splines are piecewise cubic polynomials that match first and second derivatives (as well as values)15Cubic Spline Interpolation00.10.20.30.40.50.60.70123456x valuesy valueKnown f'Nat u r a lNo KnotData16Newton Interpolating Polynomial-10123450123456X ValuesY ValuePol ynomi alData17Polynomial Applications• Data interpolation• Approximation functions in numerical quadrature and solution of ODEs• Basis functions for finite element methods• Can obtain equations for numerical differentiation• Statistical curve fitting (not discussed here) usually used in practice 18Derivative Expressions• Obtain from differentiating interpolation polynomials or from Taylor series• Series expansion for f(x) about x = a....)-(!31)-(!21)()()(333222+++−+====axdxfdaxdxfdaxdxdfafxfaxaxax∑∞===0)-(!1)(nnaxnnaxdxfdnxf• Note: d0f/dx0= f and 0! = 1• What is error from truncating series?419Truncation Error• If we truncate series after m terms∑∑∞+====+=10)-(!1)-(!1)(mnnaxnnmnnaxnnaxdxfdnaxdxfdnxfTerms used Truncation error, εm• Can write truncation error as single term at unknown location (derivation based on the theorem of the mean)1111)-(!)1(1)-(!1+=++∞+==+==∑mxmmmnnaxnnmaxdxfdmaxdxfdnξε20Derivative Expressions• Look at finite-difference grid with equal spacing: h = ∆xso xi= x0+ ih• Taylor series about x = xigives f(xi+ kh) = f[x0+ (i+k)h] = fi+kin terms of f(xi) = fi.....)(!31)(!21)()(333222++++=+===khdxfdkhdxfdkhdxdfxfkhxfiiixxxxxxiiiiixxnnnixxixxidxfdfdxfdfdxdff====== ... 22'''• Compact derivative notation21Derivative Expressions II• Combine all definitions for compact series notation.....)(!31)(!21)()(333222++++=+===khdxfdkhdxfdkhdxdfxfkhxfiiixxxxxxii• Use this formula to get expansions for various grid locations about x = xiand use results to get derivative expressions.....!3)(!2)(3'''2'''++++=+khfkhfkhfffiiiiki22Derivative Expressions III• Apply general equation for k = 1 and k = –1.....!3!23'''2'''1++++=+hfhfhfffiiiiiAhhffhfhfhfffiiiiiii+−=−−−−=++ 1'''''1'.....!3!2.!3)(!2)(3'''2'''++++=+khfkhfkhfffiiiiki.....!3!23'''2'''1+−+−=−hfhfhfffiiiiiAhhffhfhfhfffiiiiiii+−=−−+−=−−