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

1 / 55

# Bruce Mayer, PE Licensed Electrical & Mechanical Engineer BMayer@ChabotCollege.edu - PowerPoint PPT Presentation

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)

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

## PowerPoint Slideshow about 'Bruce Mayer, PE Licensed Electrical & Mechanical Engineer BMayer@ChabotCollege.edu' - kyoko

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

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 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
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 β
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 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 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:
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.:
• 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
NotationNumerical 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 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
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 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
AndAnalytical 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 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:
• 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

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 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
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 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
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

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 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
• 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 Equation2nd 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 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
MATLAB session2nd 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 Solution2nd Order ODE SUMMARY-1
• The Soln to the Homog. Eqn Produces the Complementary Solution, yc
• Assume yc take this form
• Next HOMOGENIZE the ODE
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)
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
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
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
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