W A T K I N S J O H N S O N C O M P A N Y Semiconductor Equipment Group

1
Engr/Math/Physics 25
Chp8 ODE SolnsBy MATLAB
Bruce Mayer, PE Licensed Electrical Mechanical
EngineerBMayer_at_ChabotCollege.edu
2
Learning Goals

• Use MATLABs ODE Solvers to find Solutions to
  Ordinary Differential Eqns
• When Possible use ANALYTICAL, ODE Solutions to
  Check the ACCURACY of the MATLAB Numerical Soln
• Make Approximations to perform a REALITY CHEC on
  MATLAB Solutions to NONLinear ODEs

3
MATLAB ODE Solvers

• MATLAB has several 1st Order ODE Solver Commands
  
• ode23
• ode45 (Best OverAll)
• ode113
• These Solvers based on the Runge-Kutta method,
  which is usually the best technique to apply as a
  first try for most problems.

4
MATLAB ODE Format - 1

• MATLAB ODE Solver Form

5
MATLAB ODE Format - 1

• In Vector Form (saves Writing time does NOT make
  Solution Easier)

• If we have written ROW vectors for the y b
  quantities

6
MATLAB Solution

• To solve this system using MATLAB the functions
  f1, f2, , fm must be provided to the computer,
  along with the initial values of the variables
  i.e., t0 and b1,b2, , bm
• The functions f1, f2, , fn are input using a
  function which we have to name, say fcn_vec. Then
  stored in an m-file called in this case fcn_vec.m
• The initial values of the ys are stored in a one
  dimensional array, say y0

7
MATLAB Syntax

• Obtain the Solution by typing in the command
  window

t,yode45(fcn_vec,Trng,y0)

• Where Trng ? is a vector containing at least the
  initial and final values for t.
• Example Trng T0 Tf
• Where T0 set to be t0 and Tf is set to be the
  final time at the end of the interval of
  interest.

8
Plotting the Solution

• MATLAB can also be used to plot the ODE Soluntion
  results.
• For example plot(t,y) gives a plot of ALL
  components of the solution y1, y2, , ym , as a
  function of t
• Alternatively plot(t,y(,1))gives a plot of y1 as
  a function of t

9
Comments on ode45

• Note that in its simplest form ode45 chooses its
  OWN time step and VARIES the time step according
  to how fast the solution is changing.
• Thus ode45 generates solution values at a
  sequence of times t1, t2, t3, given by tk1
  tk?tk, with ?tk selected by ode45.
• Thus EACH component of the solution y1, y2, , yn
  is ITSELF a vector containing values such as
  y1(t1), y1(t2), y1(t3), , then y2(t1), y2(t2),
  y2(t3), , etc.

10
Example ode45 (1)

• Solve 2nd ORDER ODE

• Note that to solve a 2nd order eqn we need to
  know the SLOPE (dy/dt) at t 0
• To use the 1st Order Solver, cast this 2nd Order
  eqn into 1st order (state Var) form

• With ICs

11
Example ode45 (2)

• With the Xform

• Subbing in the x1 x2

• Thus the 1st Order Eqn System in 2 Vars

• Find that

• ReArranging the ODE

12
Example ode45 (3)

• The Function File xdot_lec23.m

function xd xdot_lec23(t_val,y_vals) Bruce
Mayer, PE 05Nov05 ENGR25 Lec23 on MATLAB
ODE solvers This is the function that makes
up the system of differential equations solved
by ode45 the Vector y_vals contains yk
dy/dtk   Note the next line is necessary
because,though MATLAB does not need variables
to be initialised before they are used, vectors
will automatically be set as ROW vectors. We
want yd in COLUMN Form xdzeros(length(y_vals),1)
DEBUG Section Ttest t_val y_vals_test
y_vals xd   xd(1)y_vals(2) at t0, xdot(1)
dy(0)/dt xd(2)-5y_vals(1) -
2.0y_vals(2) sin(t_val)
13
Example ode45 (4)

• The Calling Script file Demo_ODE_Lec23.m

Bruce Mayer, PE 05Nov05 ENGR25 Lec23 on
MATLAB ODE solvers This script file calls
FUNCTION xdot_lec23 clear clear
memory set the IC's y(0) dy(0)/dt as COL
Vector y00 0 use a Time Interval of 20
Time-Units trng 0, 20   Call the ode45
routine with the above data inputs t,yode45('xd
ot_lec23', trng, y0)   Plot the first column
of the solution matrix giving x1 or
y. plot(t,y(,1)), xlabel('t'), ylabel('ODE
Solution y(t)'),... title('ODE Example -
Lecture23'), grid
14
ODE Example Result
15
Demo Problem 8.30

• Accelerating Pendulum

• For an Arbitrary Lateral-Acceleration Function,
  a(t), the ANGULAR Position, ?, described by the
  (nastily) NONlinear 2nd Order, Homogeneous ODE

q
m
mg
Time For Live Demo

• Solve for ?(t)

16
All Done for Today
FoucaultPendulum
While our clocks are set by an average 24 hour
day for the passage of the Sun from noon to noon,
the Earth rotates on its axis in 23 hours 56
minutes and 4.1 seconds with respect to the rest
of the universe. From our perspective here on
Earth, it appears that the entire universe
circles us in this time. It is possible to do
some rather simple experiments that demonstrate
that it is really the rotation of the Earth that
makes this daily motion occur.
In 1851 Leon Foucault (1819-1868) was made famous
when he devised an experiment with a pendulum
that demonstrated the rotation of the Earth..
Inside the dome of the Pantheon of Paris he
suspended an iron ball about 1 foot in diameter
from a wire more than 200 feet long. The ball
could easily swing back and forth more than 12
feet. Just under it he built a circular ring on
which he placed a ridge of sand. A pin attached
to the ball would scrape sand away each time the
ball passed by. The ball was drawn to the side
and held in place by a cord until it was
absolutely still. The cord was burned to start
the pendulum swinging in a perfect plane. Swing
after swing the plane of the pendulum turned
slowly because the floor of the Pantheon was
moving under the pendulum.
17
Engr/Math/Physics 25
Appendix
Time For Live Demo
Bruce Mayer, PE Licensed Electrical Mechanical
EngineerBMayer_at_ChabotCollege.edu
18
Prob 8.30 Results for Case-a
19
Prob 8.30 Results for Case-b
20
Prob 8.30 Results for Case-c
21
Prob 8.30 Script File
Bruce Mayer, PE 05Nov05 ENGR25 problen
8.30 This script file calls FUNCTION
pendacc clear clears memory global m b
globalize accel calc constants tmax 10 solve
over 0-20 sec tine frame set the case consts,
and IC's y(0) dy(0)/dt note that a(t) mt
b m 0, b 5, y0 0.5 0 case-a m 0,
b 5, y0 3 0 case-b m 0.5, b 0, y0
3 0 case-c Call the ode45 routine with
the above data inputs t,xode45('pendacc', 0,
tmax, y0)   Plot the first column of the
solution matrix giving x1 or y ( theta in
this case) plot(t,x(,1)), xlabel('t (sec)'),
ylabel('theta (rads)'),... title('P8.30 -
Accelerating Pendulum'), grid
22
Prob 8.30 Function File
function dxdt pendacc(t_val,z) Bruce Mayer,
PE 05Nov05 ENGR25 Prob-30 This is the
function that makes up the system of
differential equations solved by ode45 the
Vector z contains yk dy/dtk Globalize the
Constants uses to calc the Accel global m b
set the physical constants L 1 in m g
9.81 in m/sq-Sec DEBUG t_val, z  
Calc the Cauchy (State) values dxdt(1) z(2)
at t0, dxdt(1) dy(0)/dt dxdt(2)
((mt_val b)cos(z(1)) - gsin(z(1)))/L
make the dxdy into a COLUMN vector dxdt
dxdt(1) dxdt(2)