550 likes | 656 Views
Engr/Math/Physics 25. Chp9: ODE’s Numerical Solns. Bruce Mayer, PE Licensed Electrical & Mechanical Engineer BMayer@ChabotCollege.edu. Learning Goals. List Characteristics of Linear, MultiOrder, NonHomgeneous Ordinary Differential Equations (ODEs)
E N D
Engr/Math/Physics 25 Chp9: ODE’sNumerical Solns Bruce Mayer, PE Licensed Electrical & Mechanical EngineerBMayer@ChabotCollege.edu
Learning Goals • List Characteristics of Linear, MultiOrder, NonHomgeneous Ordinary Differential Equations (ODEs) • Understand the “Finite-Difference” concept that is Basis for All Numerical ODE Solvers • Use MATLAB to determine Numerical Solutions to Ordinary Differential Equations (ODEs)
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 • NONhomogenous→ RHS 0; • i.e., y(t) has a FORCING Fcn f(t) • has CONSTANT CoEfficients
Given the Simple ODE with No Zero Order (i.e., “y”) term An INITIAL Condition Solving 1st Order ODEs - 1 • AND Integrating Both Sides • Now use the IC in the Limits of Integ. • Can Solve by SEPARATING the VARIABLES • Note the use of DUMMY VARIABLES of INTEGRATIONα and β
Integrating Solving 1st Order ODEs - 2 • Separating The Variables sometimes works for 1st Order Eqns • The Function on the RHS of the 1st Order ODE is the FORCING Function • Function only of t • Can be a CONSTANT
Consider the 1st Order ODE with a “Zero” Order Term and a Forcing Fcn Solving 1st Order ODEs - 3 • yp(t) ANY Solution to the General ODE • Called the “Particular” Solution • yc(t) The Solution to the General Eqn with f(t) = 0 • The “Complementary Solution” or the “Natural” (UnForced) Response • i.e., yc is the Soln to the “Homogenous” Eqn • This is the GENERAL Eqn • By Theorems of Linear ODEs Let
Given yp and yc then the TOTAL Solution to the ODE 1st Order Response Eqns • Consider the Case Where the Forcing Function is a Constant • f(t) = A • Now Solve the ODE in Two Parts for yp & yc • For the Particular Soln, Notice that a CONSTANT Fits the Eqn:
Sub Into the General (Particular) Eqn yp and dyp/dt 1st Order Response Eqns cont • Next Separate the Variables & Integrate • Recognize LHS as a Natural Log; so • Next, Divide the Homogeneous Eqn by ·yc to yield (on whtbd) • Next Take “e” to The Power of the LHS & RHS
Then 1st Order Response Eqns cont • For This Solution Examine Extreme Cases • t = 0 • t → • is called the TIME CONSTANT • Thus the Solution for a Constant Forcing Fcn • The Latter Case is Called the Steady-State Response • All Time-Dependent Behavior has dissipated
Higher Order, Linear ODE’s • The GENERAL Higher Order ODE • Where the derivative CoEfficients, the gi(t), may be constants, including Zero • IF an analytical Solution Exists Then use the same “linear” methodology as for the First Order Eqn
Higher Order, Linear ODE’s • Where as Before • yc(t) is the solution to Complementary Eqn • yp(t) is ANY single solution to the FULL, OrginalEqn that includes the Force-Fcn • e.g.:
For More Info On Higher Order • Hi-Order ODEs usually do NOT have Analytical solns, except in special cases • Consider a 2nd order, Linear, NonHomogenous, Constant CoEfficient ODE of the form • ODE’s with these SPECIFIC Characteristics can ALWAYS be Solved Analytically • See APPENDIX for more details • These Methods used in ENGR43
Numerical ODE Solutions • Today we’ll do some MTH25 • We’ll “look under the hood” of NUMERICAL Solutions to ODE’s • The BASIC Game-Plan for even the most Sophisticated Solvers: • Given a STARTING POINT, y(0) • Use ODE to find the slopedy/dtat t=0 • ESTIMATE y1 as
Notation Numerical Solution - 1 • Exact Numerical Method (impossible to achieve) by Forward Steps yn+1 yn • Now Consider slope t tn tn+1 Dt
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
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
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
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
Euler Exmpl Calc Slope Plot
Euler vs Analytical • The Analytical Solution
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
And Analytical Soln • Now use IC • The Analytical Soln • Thus Soln u(t) • Sub u = 1−y
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 BOTH sides of the Eqn...
Use Two Steps to estimate yn+1 First → PREDICT* Use standard Euler Method to Predict Predictor-Corrector - 2 • Then Correct by using y* in the AvgCalc • Then Start the “Forward March” with the Initial Conditions
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
Predictor-Corrector Example Slope Slope
Predictor-Corrector • Greatly Improved Accuracy
ODE Example: • Euler Solution with ∆t = 0.25, y(t=0) = 37 • The Solution Table
Compare Euler vs. ODE45 Euler Solution Euler is Much LESS accurate ODE45 Solution
Compare Again with ∆t = 0.025 Euler Solution Smaller ∆T greatly improves Result ODE45 Solution
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)')
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')
All Done for Today CarlRunge Carl David Tolmé Runge Born: 1856 in Bremen, Germany Died: 1927 in Göttingen, Germany
Engr/Math/Physics 25 Appendix Time For Live Demo Bruce Mayer, PE Licensed Electrical & Mechanical EngineerBMayer@ChabotCollege.edu
Need Solutions to the 2nd Order ODE 2nd Order Linear Equation • If the Forcing Fcn is a Constant, A, Then Discern a Particular Soln • As Before The Solution Should Take This form • Verify yp • Where • yp Particular Solution • yc Complementary Solution • For Any const Forcing Fcn, f(t) = A
The Complementary Solution Satisfies the HOMOGENOUS Eqn The Complementary Solution • Look for Solution of this type • Sub Assumed Solution (y = Gest) into the Homogenous Eqn • Need yc So That the “0th”, 1st & 2nd Derivatives Have the SAME FORM so they will CANCEL (i.e., Divide-Out) in the Homogeneous Eqn • Canceling Gest • The Above is Called the Characteristic Equation
A value for “s” That SATISFIES the CHARACTERISTIC Eqn ensures that Gest is a SOLUTION to the Homogeneous Eqn Recall the Homogeneous Eqn Complementary Solution cont • The Characteristic Eqn • Solve For s by Quadratic Eqn • In terms of the Discriminant γ • If Gest is indeed a Solution Then Need
Given the “Roots” of the Homogeneous Eqn Complementary Solution cont.2 • In the Unstable case the response will grow exponentially toward ∞ • This is not terribly interesting • If the Solution is Stable, need to Consider three Sets of values for s based on the sign of γ • 1. γ > 0 → s1, s2 REAL and UNequal roots • 2. γ = 0 → s1 = s2 = s; ONE REAL root • 3. γ < 0 → Two roots as COMPLEX CONJUGATES • Can Generate STABLE and UNstable Responses • Stable • UNstable
For the Linear, 2nd Order, Constant Coeff, Homogenous Eqn Complementary Soln Cases 1&2 • By the Methods of MTH4 & ENGR43 Find Solutions to the ODE by discriminant case: Real & Unequal Roots (Stable for Neg Roots) Single Real Root (Stable for Neg Root)
Complementary Soln Case - 3 Complex Conjugate Roots of the form: s = a ± jω (Stable for Neg a) • Using the Euler Identity:And Collecting Terms find • a, ω, B1, B2 all Constants(a & ω are KNOWN)
For the Linear, 2nd Order, Constant Coeff, Homogenous Eqn 2nd Order Solution • To Find the Values of the Constants Need TWO Initial Conditions (ICs) • The ZERO Order IC • Can Find Solution based Upon the nature of the Roots of the Characteristic Eqn • The 1st Order IC
Properly Apply Initial conditions • The IC’s Apply ONLY to the TOTAL Solution • Many times It’s EASY to forget to add the PARTICULAR solution BEFORE applying the IC’s • Do NOT neglect yp(t) prior to IC’s
The Homogeneous Equation 2nd Order ODE Example - 1 • The Characteristic Eqn and Roots • And the IC’s • Then the Soln Form Given Real & UnEqual Roots
From the Zero Order IC 2nd Order ODE Example - 2 • Then at t = 0 • Now Have 2 Eqns for A1 & A2 • To Use the 1st Order IC need to take Derivative • Solve w/ MATLABBackDivision
MATLAB session 2nd Order ODE Example - 3 • The Response Curve >> C = [1,1; -3,-6]; >> b = [0.5; -8.5]; >> A = C\b A = -1.8333 2.3333 >> A_6 = 6*A A_6 = -11.0000 14.0000 Be sure to check for correct IC’s Starting-Value & Slope • Or
2nd Order ODE SuperSUMMARY-1 • See Appendix for FULL Summary • Find ANY Particular Solution to the ODE, yp (often a CONSTANT) • Homogenize ODE → set RHS = 0 • Assume yc = Gest; Sub into ODE • Find Characteristic Eqn for yc; a 2nd order Polynomial
2nd Order ODE SuperSUMMARY-2 • Find Roots to Char Eqn Using Quadratic Formula (or MATLAB) • Examine Nature of Roots to Reveal form of the Eqn for the Complementary Solution: • Real & Unequal Roots → yc = Decaying Constants • Real & Equal Roots → yc = Decaying Line • Complex Roots → yc = Decaying Sinusoid
2nd Order ODE SuperSUMMARY-3 • Then the TOTAL Solution: y = yc + yp • All TOTAL Solutions for y(t) include 2 Unknown Constants • Use the Two INITIAL Conditions to generate two Eqns for the 2 unknowns • Solve the TotalSolution for the 2 Unknowns to Complete the Solution Process
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