Bruce Mayer, PE Licensed Electrical & Mechanical Engineer BMayer@ChabotCollege.edu

1. Engr/Math/Physics 25 Chp9: ODE’sNumerical Solns Bruce Mayer, PE Licensed Electrical & Mechanical EngineerBMayer@ChabotCollege.edu

2. Learning Goals • List Characteristics of Linear, MultiOrder, NonHomgeneous Ordinary Differential Equations (ODEs) • Solve ANALYTICALLLY, Linear, 2nd Order, NonHomogeneous, Constant Coefficient ODEs • Use MATLAB to determine Numerical Solutions to Ordinary Differential Equations (ODEs)

3. Ordinary Diff Eqn Differential Equations • Partial Diff Eqn • PDE’s Not Covered in ENGR25 • Discussed in More Detail in ENGR45 • Examining the ODE, Note that it: • is LINEAR → y, dy/dt, d2y/dt2 all raised to Power of 1 • 2nd ORDER → Highest Derivative is 2 • is NONhomogenous → RHS  0; • i.e., y(t) has a FORCING Fcn f(t) • has CONSTANT Coefficients

4. Numerical ODE Solution • Today we’ll do some MTH25 • We’ll “look under the hood” of NUMERICAL Solutions to ODE’s • The BASIC Game-Plan for even Sophisticated Solvers: • Given a STARTING POINT, y(0) • Use ODE to find dy/dt at t=0 • ESTIMATE y1 as

5. Notation Numerical Solution - 1 • Exact Numerical Method (impossible to achieve) by Forward Steps yn+1 yn • Now Consider t tn tn+1 Dt

6. Numerical Solution - 2 • The diagram at Left shows that the relationship between yn, yn+1 and the CHORD slope yn+1 Tangent Slope yn Chord Slope • The problem with this formula is we canNOT calculate the chord slope exactly • We Know Only Δt & yn, but NOT the NEXT Step yn+1 t tn tn+1 Dt The AnalystChooses Δt

7. However, we can calculate the TANGENT slope at any point FROM the differential equation itself Numerical Solution -3 • The Basic Concept for all numerical methods for solving ODE’s is to use the TANGENT slope, available from the R.H.S. of the ODE, to approximate the chord slope • Recognize dy/dt as the Tangent Slope

8. Solve 1st Order ODE with I.C. Euler Method – 1st Order • ReArranging • Use: [Chord Slope]  [Tangent Slope at start of time step] • Then Start the “Forward March” with Initial Conditions

9. Consider 1st Order ODE with I.C. Euler Example • But from ODE • So In This Example: • Use The Euler Forward-Step Reln • See Next Slide for the 1st Nine Steps For Δt = 0.1

10. Euler Exmpl Calc Slope Plot

11. Euler vs Analytical • The Analytical Solution

12. Let u = −y+1 Then Analytical Soln • Integrate Both Sides • Recognize LHS as Natural Log • Sub for y & dy in ODE • Raise “e” to the power of both sides • Separate Variables

13. And Analytical Soln • Now use IC • The Analytical Soln • Thus Soln u(t) • Sub u = 1−y

14. Again Solve 1st Order ODE with I.C. Predictor-Corrector - 1 • Mathematically Avg of the Tangent Slopes at (tn,yn) & (tn+1,yn+1) • This Time Let: Chord slope average of tangent slopes at start and END of time step • BUT, we do NOT know yn+1 and it appears on the RHS...

15. Use Two Steps to estimate yn+1 First → PREDICT* Predictor-Corrector - 2 • Use y* in the Avg Calc • Then Start the “Forward March” with the Initial Conditions • Then Correct

16. Solve ODE with IC Predictor-Corrector Example • The Corrector step • The next Step Eqn for dy/dt = f(t,y)= –y+1 • Numerical Results on Next Slide

17. Predictor-Corrector Example Slope Slope

18. Predictor-Corrector • Greatly Improved Accuracy

19. ODE Example: • Euler Solution with ∆t = 0.25 • The Solution Table

20. Compare Euler vs. ODE45 Euler Solution Euler is Much LESS accurate ODE45 Solution

21. Compare Again with ∆t = 0.025 Euler Solution Smaller ∆T greatly improves Result ODE45 Solution

22. MatLAB Code for Euler % Bruce Mayer, PE % ENGR25 * 04Jan11 % file = Euler_ODE_Numerical_Example_1201.m % y0= 37; delt = 0.25; t= [0:delt:10]; n = length(t); yp(1) = y0; % vector/array indices MUST start at 1 tp(1) = 0; for k = 1:(n-1) % fence-post adjustment to start at 0 dydt = 3.9*cos(4.2*yp(k))^2-log(5.1*tp(k)+6); dydtp(k) = dydt % keep track of tangent slope tp(k+1) = tp(k) + delt; dely = delt*dydt delyp(k) = dely yp(k+1) = yp(k) + dely; end plot(tp,yp, 'LineWidth', 3), grid, xlabel('t'),ylabel('y(t) by Euler'),... title('Euler Solution to dy/dt = 3.9cos(4.2y)-ln(5.1t+6)')

23. MatLAB Command Window forODE45 >> dydtfcn = @(tf,yf) 3.9*(cos(4.2*yf))^2-log(5.1*tf+6); >> [T,Y] = ode45(dydtfcn,[0 10],[37]); >> plot(T,Y, 'LineWidth', 3), grid, xlabel('T by ODE45'), ylabel('Y by ODE45')

24. All Done for Today CarlRunge Carl David Tolmé Runge Born: 1856 in Bremen, Germany Died: 1927 in Göttingen, Germany

25. Engr/Math/Physics 25 Appendix Time For Live Demo Bruce Mayer, PE Licensed Electrical & Mechanical EngineerBMayer@ChabotCollege.edu

26. If NonHomogeneous Then find ANY Particular Solution 2nd Order ODE SUMMARY-1 • The Soln to the Homog. Eqn Produces the Complementary Solution, yc • Assume yc take this form • Next HOMOGENIZE the ODE

27. Subbing yc = Aest into the Homog. Eqn yields the Characteristic Eqn 2nd Order ODE SUMMARY-2 • Check FORM of Roots • If s1 & s2→ REAL & UNequal • Find the TWO roots that satisfy the Char Eqn by Quadratic Formula • Decaying Contant(s)

28. If s1 & s2→ REAL & Equal, then s1 = s2 =s 2nd Order ODE SUMMARY-3 • Decaying Sinusoid • Add Particlular & Complementary Solutions to yield the Complete Solution • Decaying Line • If s1 & s2→ Complex Conjugates then

29. To Find Constant Sets: (G1, G2), (m, b), (B1, B2) Take for COMPLETE solution 2nd Order ODE SUMMARY-4 • Find Number-Values for the constants to complete the solution process • Yields 2 eqns in 2 for the 2 Unknown Constants

30. Another way of thinking about numerical methods is in terms of finite differences. Use the Approximation Finite Difference Methods - 1 • And From the Differential Eqn • From these two equations obtain: • Recognize as the Euler Method

31. Could make More Accurate by Approximating dy/dt at the Half-Step as the average of the end pts Finite Difference Methods - 2 • Then Again Use the ODE to Obtain • Recognize as the Predictor-Corrector Method