Introduction to Solving Differential Equations using Matlab

1
Introduction to Solving Differential Equations
using Matlab

• Objective
• Demonstrate how the built-in Matlab function
  ode45 provides solutions to ordinary
  differential equations

2
What is a Differential Equation?

• A differential equation is an equation that
  includes derivatives of the unknown function.
• If the derivatives are ordinary derivatives (as
  opposed to partial derivatives), the equation is
  called an ordinary differential equation (ODE).
• Matlab can be used to solve several ODEs
  simultaneously
• In our examples, we will consider functions of
  time.
• Solution of the differential equation requires
  that the initial conditions (values of the
  unknown functions at time 0) be specified.

3
Examples

• Solution to a single ODE (straight line motion)
• Solution to a system of three ODEs (2D
  trajectory)
• Solution to a system of four more complicated
  ODEs (orbital mechanics)
• Pendulum motion (1st team assignment)

4
Example 1 A single ODE
Suppose that you as you drive down a highway at
55 mile/hr, you pass mile marker 100 at 1200
noon. Assuming that mile markers are increasing
in your direction of travel and you continue to
travel at 55 mph, what mile marker will you pass
at 100 PM?
5
Example 1 Analytical Solution
x(t) position (in miles) along highway as a
function of time (in hours) Governing
differential equation Initial condition x(0)
100 Solution x(t) 55 t 100 x(1) 155
6
Example 1 Matlab Solution
ode_ex1.m
ode_ex1_equations.m
7
Matlab ode45 function
time,fode45(_at_ode_ex1_equations,t,f)
output values time- array of time values f -
array of solution values
function that provides fp (df/dts)
input values t - array of times where solution
is desired f - initial value for f (unknown
variable(s))
8
Example 1 Matlab Output
Command Window
9
Example 2 System of three ODEs
Determine the trajectory of a projectile launched
at an initial velocity, v0, and initial angle, ?
10
Example 2 - Analytical Solution
Definitions t - time, x(t) - horizontal
distance, y(t) - vertical distance, g 9.81
m/s2 (accel. due to gravity) Governing
differential equations dx/dt v0 cos ? d2y/dt2
-g Initial conditions x(0) y(0) 0 v0x v0
cos ?, v0y v0 sin ? Analytical Solution x(t)
v0x t y(t) v0y t ½ g t2 tmax 2 v0y / g
11
Formulation for Matlab Solution
Define f1x(t) f2y(t) f3dy/dt Taking
derivatives with respect to time fp1 df1/dt
dx/dt v0x fp2 df2/dt dy/dt f3 fp3
df3/dt d2y/dt2 -g
12
Example 2 Matlab Solution
ode_ex2.m
13
Example 2 Matlab Solution(cont.)
ode_ex2_equations.m
14
Example 2 Matlab Output
Command Window
15
Example 3 System of four ODEs

• Determine the trajectory of a body as it orbits
  around a central body. Given the closest
  distance between the bodies (perihelion) and the
  velocity at that point, compute and animate the
  trajectory.
• Examples
• - a satellite orbiting the earth
• a comet orbiting the sun

16
Example 3 Definitions
r(t)
vo
elliptical trajectory
?(t)
R
ro
17
Example 3 Definitions (cont.)
Newtons law of gravitational attraction F G M
m / r2 where G 66.7e-12 m3/kgs2 M - mass of
central body m - mass of orbiter Non-dimensional
distance r / R Non-dimensional time T t /
?(R3 / G M)
18
Example 3 Governing ODEs
m d2r/dt2-r (d?/dt)2 -G M m / r2 m (r d2?/dt2
2 dr/dt d?/dt) 0
Initial Conditions at time0 r r0 dr/dt 0 ?
0 d?/dt v0 / r0
19
Example 3 Resulting differential equations
df1/dT f2 df2 /dT f1 f42 - 1/f12 df3 /dT
f4 df4 /dT -2 f2 f4 / f1
f1 r / R f2 (1/R) dr/dt f3 ? f4 d?/dT
gt
Initial Conditions
f1(0) ro / R f2(0) 0 f3(0) 0 f4(0) T v0 /
ro
20
Features of orbital_mechanics.m
Data input files case1.dat (satellite orbiting
the earth) case2.dat (Halleys comet orbiting
the sun)
21
Features of orbital_mechanics.m(cont.)
Main program orbital_mechanics.m User
functions orbital_solve_ode.m - solves
differential equation orbital_equations.m -
provides fp values for ode45 orbital_animation.m
- animates motion
22
Pendulum MotionTeam Programming Activity 1 -
due 4/24
Solve and animate the motion of a simple pendulum
Governing equations d2?/dt2 -(g/l) sin
? Initial conditions at time0 ?(0) ?0 d?/dt
0
23
Pendulum Motion (cont.)
Initial Conditions
f1(0) ?o f2(0) 0
24
Pendulum Motion (cont.)
Develop a program that 1. Requests user to input
initial angle, ?o 2. Generate a plot of ? vs.
time, comparing the numerical solution to the
analytical solution given by ?(t) ?o cos
sqrt(g / l) t, where l 1 m and g
9.81m/s2 3. Demonstrate your program for the
cases ?o ? / 10, ? / 4, and ? / 2