MATLAB Ordinary Differential Equations – Part II - PowerPoint Presentation

Slide1

MATLABOrdinary Differential Equations – Part II

Greg Reese,

Ph.D

Research Computing Support Group

Academic Technology Services

Miami University

September 2013Slide2

MATLABOrdinary Differential Equations – Part II

© 2010-2013 Greg Reese. All rights reserved

2Slide3

Parametric curves

3

Usually describe curve in plane with equation in

x and y, e.g.,

x

2

+ y

2

= r

2

is circle at origin

In other words, can write

y

as a function of

x

or vice-versa.

Problems

Form may not be easy to work with

Lots of curves that can't write as single equation in only

x

and

ySlide4

Parametric curves

4

One solution is to use a

parametric equation. In it, define both x and

y

to be functions of a third variable, say

t

x = f(t) y = g(t)

Each value of

t

defines a point

(

x,y

)=( f(t), g(t) )

that we can plot. Collection of points we get by letting

t

take on all its values is a

parametric curveSlide5

Parametric curves

5

To plot 2D parametric curves use

ezplot

(

funx,funy

)

ezplot

(

funx,funy

,[

tmin,tmax

])

where

funx

is handle to function

x(t)

funy

is handle to function

y(t)

tmin

,

tmax

specify range of

t

if omitted, range is

0 < t < 2

πSlide6

Parametric curves

6

Try It

Plot x(t) = cos

(t) y(t) = 3sin(t)

over

0<t<2

π

>>

ezplot

( @(t)

cos

(t), @(t)3*sin(t) )Slide7

Parametric curves

7

Try It

Plot x(t) = cos

3

(t) y(t) = 3sin(t)

Hint: use the vector arithmetic operators

.*, ./, and .^ to avoid warnings

>>

ezplot

( @(t)

cos

(t).^3, @(t)3*sin(t) )Slide8

Parametric curves

8

Often parametric curves expressed in polar form

ρ = f(θ)

Plot with

ezpolar

(fun)

ezpolar

(f,[

thetaMin,thetaMax

])

where

f

is handle to function

ρ

= f(

θ

)

thetaMin

,

thetaMax

specify range of

θ

if omitted, range is

0 <

θ < 2π

θ

ρ

(

θ

)

x

ySlide9

Parametric curves

9

Try It

Plot ρ =

θ

>>

ezpolar

( @(theta)theta )Slide10

Parametric curves

10

Try It

Plot ρ

= sin(

θ

)

cos

(

θ

)

over

0 <

θ

< 6

π

Hint: use the vector arithmetic operators .*, ./, and .^ to avoid warnings

>> ezpolar( @(theta)sin(theta).*cos(theta),

...

[0 6*pi ] )Slide11

Parametric curves

11

To plot 3D parametric curves use

ezplot

(

funx,funy,funz

)

ezplot

(

funx,funy,funz

,[

tmin,tmax

])

where

funx

is handle to function

x(t)

funy

is handle to function

y(t)

funz

is handle to function

z(t)

tmin,

tmax specify range of

t if omitted, range is 0 < t < 2πSlide12

Parametric curves

12

Try It

Plot

x(t) =

cos

(4t) y(t) = sin(4t) z(t) = t

over

0<t<2

π

>> ezplot3( @(t)

cos

(4*t), @(t)sin(4*t), @(t)t )

Unfortunately, there's no ezpolar3Slide13

Parametric curves

Try It

Plot

x(t) =

cos

(4t) y(t) = sin(4t) z(t) = t

over

0<t<2

π

>>ezplot3( @(t)

cos

(4*t),@(t)sin(4*t),... @(t)t ), 'animate ' )

13

FOR THRILLSSlide14

Try It

Place command window and figure window side by side and use

comet3()

to plot x(t) = cos(30t) y(t) = sin(30t) z(t) = t over 0<t<2π

>> t = 2*pi * 0:0.001:1;

>> x =

cos

( 30*t );

>> y = sin( 30*t );

>> z = t;

>> comet3( x, y, z )

Parametric curves

14

FOR THRILLS

and

CHILLSSlide15

Parametric curves

15

Questions?Slide16

Phase plane plot

16

For ideal pendulum

, θ

''

+sin(

θ

(t) ) = 0

Define

y

1

(t) =

θ

(t), y

2

(t) =

θ'(t) to get

Write pendulum.m

function

yPrime

=

pendulum

( t, y )

yPrime

= [ y(2) -sin( y(1) ) ]';

θ(t)

gravitySlide17

Phase plane plot

17

3 different initial conditions

y

1

(0)= θ(0) = 1

R

=57°

y

2

(0)= θ

'

(0) = 1

R

/sec=57°/sec

θ(0)

θ

'

(0)

y

1

(0)= θ(0)

= -5

R

=-286°=74°

y

2

(0)= θ

'

(0) = 2

R

/sec=115°/sec

θ(0)

θ

'

(0)

y

1

(0)= θ(0)

= 5

R

=-74°

y

2

(0)= θ

'

(0) =

-2

R

/sec=-115°/sec

θ(0)

θ

'

(0)Slide18

Phase plane plot

18

Try It

For ideal pendulum

,

θ

''

+sin(

θ

(t) ) = 0

solve for the initial conditions

θ

(0)=1,

θ

'

(0)=1

and time = [ 0 10 ] and make a phase plane plot with y1

(t) on the horizontal axis and y2

(t)

on the vertical. Store the results in

ta

and

yaSlide19

Phase plane plot

19

Try It

>>

tSpan

= [ 0 10 ];

>> y0 = [ 1 1 ];

>> [

ta

ya

] = ...

ode45( @pendulum,

tSpan

, y0 );

>> plot(

ya

(:,1),

ya

(:,2) );

y

1

(0)= θ(0) = 1

R

=57°

y

2

(0)= θ

'

(0) = 1

R

/sec=57°/sec

θ(0)

θ

'

(0)

Qualitatively,

what should pendulum do?Slide20

Phase plane plot

20

Try It

>>

tSpan

= [ 0 10 ];

>> y0 = [ -5 2 ];

>> [

tb

yb

] = ...

ode45( @pendulum,

tSpan

, y0 );

>> plot(

yb

(:,1),

yb

(:,2) );

Qualitatively,

what should pendulum do?

y

1

(0)= θ(0)

= -5

R

=-286°=74°

y

2

(0)= θ

'

(0) = 2

R

/sec=115°/sec

θ(0)

θ

'

(0)Slide21

Phase plane plot

21

Try It

>>

tSpan

= [ 0 10 ];

>> y0 = [ 5 -2 ];

>> [

tc

yc

] = ...

ode45( @pendulum,

tSpan

, y0 );

>> plot(

yc

(:,1),

yc

(:,2) );

Qualitatively,

what should pendulum do?

y

1

(0)= θ(0)

= 5

R

=-74°

y

2

(0)= θ

'

(0) =

-2

R

/sec=-115°/sec

θ(0)

θ

'

(0)Slide22

Phase plane plot

22

Try It

Graph all three on one plot

>>

plot

( ya(:,1), ya(:,2),

yb

(:,1),

yb

(:,2),...

yc

(:,1),

yc

(:,2) )

>>

ax

= axis;

>> axis( [ -5 5

ax

(3:4) ] );Slide23

Phase plane plot

23

If have initial

condition (other

than previous 3)

that is exactly on

curve (red dot)

can tell its path in

phase plane.

Q: What if not on curve but very close to it

(yellow dot)? A: ?Slide24

Phase plane plot

24

To help understand solution for any initial condition, can make phase plot and add information about how each state variable changes with time, i.e., display the first derivative of each state variable.Slide25

Phase plane plot

25

Will show rate of change of state variables at a point by drawing a vector point there. Horizontal component of vector is rate of change of variable one; vertical component of vector is rate of change of variable two.

y

1

(t)

y

2

(t)

y

'

1

(t)

y

'

2

(t)Slide26

Phase plane plot

26

Where can we get these rates of change?

From the state-space formulation

y

'

(t) = f( t, y )

!

Example – ideal pendulum

y

1

(t)

y

2

(t)

y

'

1

(t)

y

'

2

(t)Slide27

Phase plane plot

27

To plot vectors at point, use

quiver( x, y, u, v )

This plots the vectors (

u,v

)

at every point (

x,y

)

x

is matrix of x-values of points

y

is matrix of y-values of points

u

is matrix of horizontal components of vectors

v

is matrix of vertical components of vectors

All matrices must be same size

v

(

x,y

)

uSlide28

Phase plane plot

28

To make

x and

y

for quiver, use

[ x y ] =

meshgrid

(

xVec

,

yVec

)

Example

>> [ x y ] =

meshgrid

( 1:5, 7:9 )

x = 1 2 3 4 5

1 2 3 4 5

1 2 3 4 5

y = 7 7 7 7 7

8 8 8 8 8

9 9 9 9 9Slide29

Phase plane plot

29

Let's make quiver plot at every point

(x,y

)

for

x

going from -5 to 5 in increments of 1 and

y

going from -2.5 to 2.5 in

increments of 0.5

>> [y1 y2 ] =

meshgrid

( -5:5, -2.5:0.5:2.5 );

>> y1Prime = y2;

>> y2Prime = -sin( y1 );

>>

quiver

( y1, y2, y1Prime, y2Prime )Slide30

Phase plane plot

30

Now can see rate of change of state variables. MATLAB plots zero-vectors as a small dot.

What is physical meaning of 3 small dots?Slide31

Phase plane plot

31

To answer question about solution with initial conditions close to those of another solution (yellow dot close to green line), put phase-plane plot and quiver plot togetherSlide32

Phase plane plot

32

Try It

>> plot(

ya

(:,1),

ya

(:,2),

yb

(:,1), ...

yb

(:,2),

yc

(:,1),

yc

(:,2) )

>> ax = axis;

>> axis( [ -5 5 ax(3:4) ] )

>> hold on

>> quiver( y1, y2, y1Prime, y2Prime )

>> hold offSlide33

Phase plane plot

33

Try It

To see solution path

for specific initial

conditions, imagine

dropping a toy boat (initial condition) at a spot in a river (above plot) and watching how current (arrows) pushes it around.Slide34

Phase plane plot

34

What path would dot take and why?Slide35

Phase plane plot

35

From phase-plane plot it appears reasonable to say that if the initial conditions of the solutions of a differential equation are close to each other, the solutions are also close to each other.Slide36

Phase plane plot

36

Let's check out this idea that close initial conditions lead to close solutions a little more.

Replot

solution to first initial conditions

>>

tSpan

= [ 0 10 ];

>> y0a = [ 1 1 ];

>> [

ta

ya

] = ...

ode45( @pendulum,

tSpan

, y0a );

>> plot(

ya

(:,1),

ya

(:,2) );Slide37

Phase plane plot

37

Now let's solve again with initial conditions 25% greater and plot both

>> n = 0.25;

>> yy0 = y0a + n*y0a;

>> [

tt

yy

] = ode45(@pendulum,

tSpan

, yy0 );

>> plot(

ya

(:,1),

ya

(:,2),

yy

(:,1),

yy

(:,2) )Slide38

Phase plane plot

38

Fairly closeSlide39

Phase plane plot

39

Repeat for 10% greater

>> n = 0.1;

>> yy0 = y0a + n*y0a;

>> [

tt

yy

] = ode45( @pendulum,

tSpan

, yy0 );

>> plot(

ya

(:,1),

ya

(:,2),

yy

(:,1),

yy

(:,2) )Slide40

Phase plane plot

40Slide41

Phase plane plot

41

Repeat for 1% greater

>> n = 0.01;

>> yy0 = y0a + n*y0a;

>> [

tt

yy

] = ode45( @pendulum,

tSpan

, yy0 );

>> plot(

ya

(:,1),

ya

(:,2),

yy

(:,1),

yy

(:,2) )Slide42

Phase plane plot

42Slide43

Phase plane plot

43

Repeat for 0.1% greater

>> n = 0.001;

>> yy0 = y0a + n*y0a;

>> [

tt

yy

] = ode45( @pendulum,

tSpan

, yy0 );

>> plot(

ya

(:,1),

ya

(:,2),

yy

(:,1),

yy

(:,2) )Slide44

Phase plane plot

44Slide45

Phase plane plot

45

Again, it appears that if the initial conditions of the solutions of a differential equation are close to each other, the solutions are also close to each other.Slide46

Phase plane plot

46

Well, as that famous philosopher might saySlide47

Phase plane plots

47

Questions?Slide48

48

CHAOS

or

Welcome to my WorldSlide49

Chaos

49

Chaos theory

is branch of math that studies behavior of certain kinds of dynamical systems Chaotic behavior observed in nature, e.g., weather

Quantum chaos theory

studies relationship between chaos and quantum mechanics Slide50

Chaos

50

Chaotic systems are:

Deterministic – no randomness involved

If start with

identical

initial conditions, get

identical

final states

High sensitivity to initial conditions

Tiny differences in starting state can lead to enormous differences in final state, even over small time ranges

Seemingly random

Unexpected and abrupt changes in state occur

Often sensitive to slight parameter changesSlide51

Chaos

51

In 1963, Edward Lorenz, mathematician and meteorologist, published set of equations

Simplified model of convection rolls in the atmosphereAlso used as simple model of laser and dynamo (electric generator)Slide52

Chaos

52

Set of equations

NonlinearThree-dimensionalDeterministic, i.e.,

no randomness involved

Important implications for climate and weather prediction

Atmospheres may exhibit quasi-periodic behavior and may have abrupt and

seemingly

random change, even if fully deterministic

Weather can't be predicted too far into future!Slide53

Chaos

53

Equations, in state-space form, are

*

Notice only two terms have nonlinearities

* Also appear in slightly different formsSlide54

Chaos

54

All parameters are > 0

β usually 8/3σ (Prandtl number) usually 10

ρ

(Rayleigh number) often variedSlide55

Chaos

55

Let's check out the system

Download lorenz.m

function

yPrime

= ...

lorenz

(

t,y,beta,rho,sigma

)

yPrime

= zeros( 3, 1 );

yPrime

(1) = sigma*( y(2) - y(1) );

yPrime

(2) = y(1)*( rho - y(3) ) - y(2);

yPrime

(3) = y(1)*y(2) - beta*y(3);Slide56

Chaos

56

Try

It

Look at

effect

of

changing

ρ

, start with

ρ

=12

>> beta = 8/3;

>> rho = 12;

>> sigma = 10;

>>

tSpan

= [ 0 50 ];

>> y0 = [1e-8 0 0 ];

>> [ t y ] = ode45( @(

t,y

)

lorenz

(...

t,y,beta,rho,sigma

),

tSpan

, y0 );

>> plot3( y(:,1), y(:,2), y(:,3) )

>> comet3( y(:,1), y(:,2), y(:,3) )Slide57

Chaos

57

Try

It

System

converges

Az

: -120 El: -20Slide58

Chaos

58

Try

It

View

two

state

variables at a time

>>

plot

( y(:,1), y(:,2) )

>>

plot

( y(:,1), y(:,3) )

>>

plot

( y(:,2), y(:,3) )Slide59

Chaos

59

Try

It

Set

ρ

= 16

rerun

, and do comet3, plot3

Az

: -120 El: -20Slide60

Chaos

60

Try

It

Notice that "particle" gets pulled over into "hole". "Hole" is called an

attractor

Az

: -120 El: -20Slide61

Chaos

61

Try

It

Set

ρ

= 20

rerun

, and do comet3, plot3

Az

: -120 El: -20Slide62

Chaos

62

Try

It

Set

ρ

=

24.2

rerun

, and do comet3, plot3

Az

: -120 El: -20Slide63

Chaos

63

Try

It

Set

ρ

=

24.3

rerun

, and do comet3, plot3

Watch

comet

carefully

!

Az

: -120 El: -20Slide64

Chaos

64

Try

It

Wow

! A

small

change

in

ρ

causes giant change in trajectory! Particle starts bouncing back and forth between attractors

Az

: -120 El: -20Slide65

Chaos

65

Try

It

Set

ρ

= 26

rerun

, and do comet3, plot3

Az

: -120 El: -20Slide66

Chaos

66

Try

It

Set

ρ

= 30

rerun

, and do comet3, plot3

In

comet

,

watch

hopping

back and

forth

between

attractors

Az

: -120 El: -20Slide67

Chaos

67

More

common

view

of

Lorenz

attractor

Has

been

shown

Bounded

(

always

within

a box)

Non-

periodic

Never

crosses

itself

Az

: 0 El: 0Slide68

Chaos

68

Lorenz

attractor

shows

us

some

characteristics

of

chaotic

systems

Paths

in

phase

space

can

be

very

complicated

Paths can

have

abrupt

changes at seemingly

random timesSmall

variations in a parameter can produce

large changes in trajectoriesSlide69

Chaos

69

Now

look at sensitivity

to

initial

conditions

>> beta = 8/3;

>> sigma = 10;

>>

rho = 28

;

>> y0 = 1e-8 * [ 1 0 0 ];

>> [

tt

yy

] = ode45( @(

t,y

)

lorenz

(

t,y,beta,rho,sigma

),

tSpan

, y0 );

>> plot3(

yy

(:,1),

yy(:,2),

yy(:,3), 'b' )>>

title( 'Original' )

Az

: -120 El: -20Slide70

Chaos

70

Now

look at sensitivity

to

initial

conditions

>> y =

yy

;

>> plot3(

yy

(:,1),

yy

(:,2),

yy

(:,3), 'b',... y(:,1), y(:,2), y(:,3), 'y' )

>>

title

( '0%

difference

' )

OOPS! MATLAB

bug

.

Yellow

should

exactly

cover blue…

Justpretend

it does!

Az: -120 El: -20Slide71

Chaos

71

Now

look at sensitivity

to

initial

conditions

>> y0 = 1e-8 * [

1.01

0 0 ];

>> [ t y ] = ode45( @(

t,y

)

lorenz

(

t,y,beta,rho,sigma

),

tSpan

, y0 );

>> plot3(

yy

(:,1),

yy

(:,2),

yy

(:,3), 'b', y(:,1), y(:,2), y(:,3), 'y' )

>>

title

( '1%

difference' )

(1.01-1.00)/1 x 100

= 1% differenceSlide72

Chaos

72

1%

difference –

clearly

different

paths

Az

: -120 El: -20Slide73

Chaos

73

>>y0=1e-8*[

1.001

0 0 ]; % 0.1%

difference

Az

: -120 El: -20Slide74

Chaos

74

>>y0=1e-8*[

1.00001

0 0 ]; % 0.001%

difference

Az

: -120 El: -20Slide75

Chaos

75

>>y0=1e-8*[

1.0000001

0 0 ]; % 0.00001%

difference

Az

: -120 El: -20Slide76

Chaos

76

>>y0=1e-8*[

1.0000001

0 0 ]; % 0.00001%

difference

y

1

(t) vs y

2

(t)Slide77

Chaos

77

>>y0=1e-8*[

1.0000001

0 0 ]; % 0.00001%

difference

y

1

(t) vs y

3

(t)Slide78

Chaos

78

>>y0=1e-8*[

1.0000001

0 0 ]; % 0.00001%

difference

y

2

(t) vs y

3

(t)Slide79

Chaos

79

So

even

though

initial

conditions

only

differ

by

1 / 100,000 of a

percent

, the

trajectories

become

very

different

!Slide80

Chaos

80

This

extreme sensitivity

to

initial

conditions

is

often

called

The

Butterfly

Effect

A

butterfly

flapping

its

wings

in Brazil can cause a tornado in Texas.Slide81

Chaos

81

Lorenz equations are good example of chaotic system

Deterministic High sensitivity to initial conditions

Very tiny differences in starting state can lead to substantial differences in final state

Unexpected and abrupt changes in state

Sensitive to slight parameter changesSlide82

Chaos

82

MATLAB is good for studying chaotic systems

Easy to set and change initial conditions or parametersSolving equations is fast and easy

Plotting and comparing 2D and 3D trajectories also fast and easySlide83

Chaos

83

Questions?Slide84

84

The End