slide1 n.
Skip this Video
Loading SlideShow in 5 Seconds..
Bruce Mayer, PE Licensed Electrical & Mechanical Engineer PowerPoint Presentation
Download Presentation
Bruce Mayer, PE Licensed Electrical & Mechanical Engineer

Loading in 2 Seconds...

play fullscreen
1 / 55

Bruce Mayer, PE Licensed Electrical & Mechanical Engineer - PowerPoint PPT Presentation

  • Uploaded on

Engr/Math/Physics 25. Chp9: ODE’s Numerical Solns. Bruce Mayer, PE Licensed Electrical & Mechanical Engineer Learning Goals. List Characteristics of Linear, MultiOrder, NonHomgeneous Ordinary Differential Equations (ODEs)

I am the owner, or an agent authorized to act on behalf of the owner, of the copyrighted work described.
Download Presentation

PowerPoint Slideshow about 'Bruce Mayer, PE Licensed Electrical & Mechanical Engineer' - kyoko

An Image/Link below is provided (as is) to download presentation

Download Policy: Content on the Website is provided to you AS IS for your information and personal use and may not be sold / licensed / shared on other websites without getting consent from its author.While downloading, if for some reason you are not able to download a presentation, the publisher may have deleted the file from their server.

- - - - - - - - - - - - - - - - - - - - - - - - - - E N D - - - - - - - - - - - - - - - - - - - - - - - - - -
Presentation Transcript

Engr/Math/Physics 25

Chp9: ODE’sNumerical Solns

Bruce Mayer, PE

Licensed Electrical & Mechanical

learning goals
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)
differential equations
Ordinary Diff EqnDifferential 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
solving 1 st order odes 1
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.
  • Note the use of DUMMY VARIABLES of INTEGRATIONα and β
solving 1 st order odes 2
IntegratingSolving 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
solving 1 st order odes 3
Consider the 1st Order ODE with a “Zero” Order Term and a Forcing FcnSolving 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
1 st order response eqns
Given yp and yc then the TOTAL Solution to the ODE1st 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:
1 st order response eqns cont
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
1 st order response eqns cont1
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
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 s1
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
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
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
numerical solution 1
NotationNumerical Solution - 1
  • Exact Numerical Method (impossible to achieve) by Forward Steps



  • Now Consider slope





numerical solution 2
Numerical Solution - 2
  • The diagram at Left shows that the relationship between yn, yn+1 and the CHORD slope


Tangent 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





The AnalystChooses Δt

numerical solution 3
However, we can calculate the TANGENT slope at any point FROM the differential equation itselfNumerical 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
euler method 1 st order
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
euler example
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 vs analytical
Euler vs Analytical
  • The Analytical Solution
analytical soln
Let u = −y+1


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
analytical soln1
AndAnalytical Soln
  • Now use IC
  • The Analytical Soln
  • Thus Soln u(t)
  • Sub u = 1−y
predictor corrector 1
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...
predictor corrector 2
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
predictor corrector example
Solve ODE with ICPredictor-Corrector Example
  • The Corrector step
  • The next Step Eqn for dy/dt = f(t,y)= –y+1
  • Numerical Results on Next Slide
predictor corrector
  • Greatly Improved Accuracy
ode example
ODE Example:
  • Euler Solution with ∆t = 0.25, y(t=0) = 37
  • The Solution Table
compare euler vs ode45
Compare Euler vs. ODE45

Euler Solution

Euler is Much LESS accurate

ODE45 Solution

compare again with t 0 0 25
Compare Again with ∆t = 0.025

Euler Solution

Smaller ∆T greatly improves Result

ODE45 Solution

matlab code for euler
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;


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 for ode45
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
All Done for Today


Carl David Tolmé Runge

Born: 1856 in Bremen, Germany

Died: 1927 in Göttingen, Germany


Engr/Math/Physics 25


Time For

Live Demo

Bruce Mayer, PE

Licensed Electrical & Mechanical

2 nd order linear equation
Need Solutions to the 2nd Order ODE2nd 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
The Complementary Solution Satisfies the HOMOGENOUS EqnThe 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
complementary solution cont
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
complementary solution cont 2
Given the “Roots” of the Homogeneous EqnComplementary 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
complementary soln cases 1 2
For the Linear, 2nd Order, Constant Coeff, Homogenous EqnComplementary 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
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)
2 nd order solution
For the Linear, 2nd Order, Constant Coeff, Homogenous Eqn2nd 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
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
2 nd order ode example 1
The Homogeneous Equation2nd Order ODE Example - 1
  • The Characteristic Eqn and Roots
  • And the IC’s
  • Then the Soln Form Given Real & UnEqual Roots
2 nd order ode example 2
From the Zero Order IC2nd 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
2nd order ode example 3
MATLAB session2nd Order ODE Example - 3
  • The Response Curve

>> C = [1,1; -3,-6];

>> b = [0.5; -8.5];

>> A = C\b

A =



>> A_6 = 6*A

A_6 =



Be sure to check for correct IC’s  Starting-Value & Slope

  • Or
2 nd order ode supersummary 1
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
2 nd order ode supersummary 2
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
2 nd order ode supersummary 3
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
2 nd order ode summary 1
If NonHomogeneous Then find ANY Particular Solution2nd Order ODE SUMMARY-1
  • The Soln to the Homog. Eqn Produces the Complementary Solution, yc
  • Assume yc take this form
2 nd order ode summary 2
Subbing yc = Aest into the Homog. Eqn yields the Characteristic Eqn2nd 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)
2 nd order ode summary 3
If s1 & s2→ REAL & Equal, then s1 = s2 =s2nd Order ODE SUMMARY-3
    • Decaying Sinusoid
  • Add Particlular & Complementary Solutions to yield the Complete Solution
    • Decaying Line
  • If s1 & s2→ Complex Conjugates then
2 nd order ode summary 4
To Find Constant Sets: (G1, G2), (m, b), (B1, B2) Take for COMPLETE solution2nd 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
finite difference methods 1
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
finite difference methods 2
Could make More Accurate by Approximating dy/dt at the Half-Step as the average of the end ptsFinite Difference Methods - 2
  • Then Again Use the ODE to Obtain
  • Recognize as the Predictor-Corrector Method