MATLAB Ordinary Differential Equations Part II Greg Reese

MATLAB Ordinary Differential Equations – Part II Greg Reese, Ph. D Research Computing Support Group Academic Technology Services Miami University September 2013

MATLAB Ordinary Differential Equations – Part II © 2010 -2013 Greg Reese. All rights reserved 2

Parametric curves 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 y 3

Parametric curves 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 curve 4

Parametric curves To plot 2 D 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π 5

Parametric curves Try It Plot x(t) = cos(t) y(t) = 3 sin(t) over 0<t<2π >> ezplot( @(t)cos(t), @(t)3*sin(t) ) 6

Parametric curves Try It Plot x(t) = cos 3(t) y(t) = 3 sin(t) – Hint: use the vector arithmetic operators. *, . /, and. ^ to avoid warnings >> ezplot( @(t)cos(t). ^3, @(t)3*sin(t) ) 7

Parametric curves Often parametric curves expressed in polar form ρ = f(θ) Plot with y ρ(θ) θ x ezpolar(fun) ezpolar(f, [theta. Min, theta. Max]) where • f is handle to function ρ = f(θ) • theta. Min, theta. Max specify range of θ – if omitted, range is 0 < θ < 2π 8

Parametric curves Try It Plot ρ = θ >> ezpolar( @(theta)theta ) 9

Parametric curves 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 ] ) 10

Parametric curves To plot 3 D 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π 11

Parametric curves Try It Plot x(t) = cos(4 t) y(t) = sin(4 t) z(t) = t over 0<t<2π >> ezplot 3( @(t)cos(4*t), @(t)sin(4*t), @(t)t ) Unfortunately, there's no ezpolar 3 12

FOR THRILLS Parametric curves Try It Plot x(t) = cos(4 t) y(t) = sin(4 t) z(t) = t over 0<t<2π >>ezplot 3( @(t)cos(4*t), @(t)sin(4*t), . . . @(t)t ), 'animate' ) 13

Parametric curves Try It FOR THRILLS and CHILLS Place command window and figure window side by side and use comet 3() to plot x(t) = cos(30 t) y(t) = sin(30 t) z(t) = t over 0<t<2π >> >> >> t = 2*pi * 0: 0. 001: 1; x = cos( 30*t ); y = sin( 30*t ); z = t; comet 3( x, y, z ) 14

Parametric curves Questions? 15

Phase plane plot For ideal pendulum, θ''+sin( θ(t) ) = 0 Define y 1(t) = θ(t), y 2(t) = θ'(t) to get θ(t) gravity Write pendulum. m function y. Prime = pendulum( t, y ) y. Prime = [ y(2) -sin( y(1) ) ]'; 16

Phase plane plot 3 different initial conditions θ'(0) θ(0) θ'(0) y 1(0)= θ(0) = 1 R=57° y 2(0)= θ'(0) = 1 R/sec=57°/sec y 1(0)= θ(0) = -5 R=-286°=74° y 2(0)= θ'(0) = 2 R/sec=115°/sec y 1(0)= θ(0) = 5 R=-74° y 2(0)= θ'(0) = -2 R/sec=-115°/sec 17

Phase plane plot 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 y 1(t) on the horizontal axis and y 2(t) on the vertical. Store the results in ta and ya 18

Phase plane plot Qualitatively, what should pendulum do? Try It >> t. Span = [ 0 10 ]; >> y 0 = [ 1 1 ]; >> [ ta ya ] =. . . ode 45( @pendulum, t. Span, y 0 ); >> plot( ya(: , 1), ya(: , 2) ); θ(0) θ'(0) y 1(0)= θ(0) = 1 R=57° y 2(0)= θ'(0) = 1 R/sec=57°/sec 19

Phase plane plot θ'(0) θ(0) Qualitatively, what should pendulum do? Try It >> t. Span = [ 0 10 ]; >> y 0 = [ -5 2 ]; >> [ tb yb ] =. . . ode 45( @pendulum, t. Span, y 0 ); >> plot( yb(: , 1), yb(: , 2) ); y 1(0)= θ(0) = -5 R=-286°=74° y 2(0)= θ'(0) = 2 R/sec=115°/sec 20

Phase plane plot θ'(0) θ(0) Qualitatively, what should pendulum do? Try It >> t. Span = [ 0 10 ]; >> y 0 = [ 5 -2 ]; >> [ tc yc ] =. . . ode 45( @pendulum, t. Span, y 0 ); >> plot( yc(: , 1), yc(: , 2) ); y 1(0)= θ(0) = 5 R=-74° y 2(0)= θ'(0) = -2 R/sec=-115°/sec 21

Phase plane plot 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) ] ); 22

Phase plane plot 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: ? 23

Phase plane plot 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. 24

Phase plane plot 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) 25

Phase plane plot Where can we get these rates of change? y (t) y'1(t) y'2(t) 2 From the state-space formulation y 1(t) y'(t) = f( t, y ) ! Example – ideal pendulum 26

Phase plane plot To plot vectors at point, use quiver( x, y, u, v ) u (x, y) 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 27

Phase plane plot To make x and y for quiver, use [ x y ] = meshgrid( x. Vec, y. Vec ) Example >> [ x y ] = meshgrid( 1: 5, 7: 9 ) x = 1 2 3 4 5 y = 7 8 9 7 8 9 28

Phase plane plot 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 >> >> [y 1 y 2 ] = meshgrid( -5: 5, -2. 5: 0. 5: 2. 5 ); y 1 Prime = y 2; y 2 Prime = -sin( y 1 ); quiver( y 1, y 2, y 1 Prime, y 2 Prime ) 29

Phase plane plot 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? 30

Phase plane plot 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 together 31

Phase plane plot 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( y 1, y 2, y 1 Prime, y 2 Prime ) >> hold off 32

Phase plane plot 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. 33

Phase plane plot What path would dot take and why? 34

Phase plane plot 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. 35

Phase plane plot Let's check out this idea that close initial conditions lead to close solutions a little more. Replot solution to first initial conditions >> t. Span = [ 0 10 ]; >> y 0 a = [ 1 1 ]; >> [ ta ya ] =. . . ode 45( @pendulum, t. Span, y 0 a ); >> plot( ya(: , 1), ya(: , 2) ); 36

Phase plane plot Now let's solve again with initial conditions 25% greater and plot both >> >> n = 0. 25; yy 0 = y 0 a + n*y 0 a; [ tt yy ] = ode 45(@pendulum, t. Span, yy 0 ); plot( ya(: , 1), ya(: , 2), yy(: , 1), yy(: , 2) ) 37

Phase plane plot Fairly close 38

Phase plane plot Repeat for 10% greater >> >> n = 0. 1; yy 0 = y 0 a + n*y 0 a; [ tt yy ] = ode 45( @pendulum, t. Span, yy 0 ); plot( ya(: , 1), ya(: , 2), yy(: , 1), yy(: , 2) ) 39

Phase plane plot 40

Phase plane plot Repeat for 1% greater >> >> n = 0. 01; yy 0 = y 0 a + n*y 0 a; [ tt yy ] = ode 45( @pendulum, t. Span, yy 0 ); plot( ya(: , 1), ya(: , 2), yy(: , 1), yy(: , 2) ) 41

Phase plane plot 42

Phase plane plot Repeat for 0. 1% greater >> >> n = 0. 001; yy 0 = y 0 a + n*y 0 a; [ tt yy ] = ode 45( @pendulum, t. Span, yy 0 ); plot( ya(: , 1), ya(: , 2), yy(: , 1), yy(: , 2) ) 43

Phase plane plot 44

Phase plane plot 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. 45

Phase plane plot Well, as that famous philosopher might say 46

Phase plane plots Questions? 47

CHAOS or Welcome to my World 48

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 49

Chaos 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 changes 50

Chaos In 1963, Edward Lorenz, mathematician and meteorologist, published set of equations • Simplified model of convection rolls in the atmosphere • Also used as simple model of laser and dynamo (electric generator) 51

Chaos Set of equations • Nonlinear • Three-dimensional • Deterministic, 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! 52

Chaos Equations, in state-space form, are* Notice only two terms have nonlinearities * Also appear in slightly different forms 53

Chaos • All parameters are > 0 • β usually 8/3 • σ (Prandtl number) usually 10 • ρ (Rayleigh number) often varied 54

Chaos Let's check out the system Download lorenz. m function y. Prime =. . . lorenz(t, y, beta, rho, sigma) y. Prime = zeros( 3, 1 ); y. Prime(1) = sigma*( y(2) - y(1) ); y. Prime(2) = y(1)*( rho - y(3) ) - y(2); y. Prime(3) = y(1)*y(2) - beta*y(3); 55

Chaos Try It Look at effect of changing ρ, start with ρ=12 >> beta = 8/3; >> rho = 12; >> sigma = 10; >> t. Span = [ 0 50 ]; >> y 0 = [1 e-8 0 0 ]; >> [ t y ] = ode 45( @(t, y)lorenz(. . . t, y, beta, rho, sigma), t. Span, y 0 ); >> plot 3( y(: , 1), y(: , 2), y(: , 3) ) >> comet 3( y(: , 1), y(: , 2), y(: , 3) ) 56

Chaos Try It System converges Az: -120 El: -20 57

Chaos Try It View two state variables at a time >> plot( y(: , 1), y(: , 2) ) >> plot( y(: , 1), y(: , 3) ) >> plot( y(: , 2), y(: , 3) ) 58

Chaos Try It Set ρ = 16 rerun, and do comet 3, plot 3 Az: -120 El: -20 59

Chaos Try It Notice that "particle" gets pulled over into "hole". "Hole" is called an attractor Az: -120 El: -20 60

Chaos Try It Set ρ = 20 rerun, and do comet 3, plot 3 Az: -120 El: -20 61

Chaos Try It Set ρ = 24. 2 rerun, and do comet 3, plot 3 Az: -120 El: -20 62

Chaos Try It Set ρ = 24. 3 rerun, and do comet 3, plot 3 Watch comet carefully! Az: -120 El: -20 63

Chaos Try It Wow! A small change in ρ causes giant change in trajectory! Particle starts bouncing back and forth between attractors 64 Az: -120 El: -20

Chaos Try It Set ρ = 26 rerun, and do comet 3, plot 3 Az: -120 El: -20 65

Try It Chaos Set ρ = 30 rerun, and do comet 3, plot 3 In comet, watch hopping back and forth between attractors Az: -120 El: -20 66

Chaos More common view of Lorenz attractor Has been shown – Bounded (always within a box) – Non-periodic – Never crosses itself Az: 0 El: 0 67

Chaos 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 times • Small variations in a parameter can produce large changes in trajectories 68

Chaos Now look at sensitivity to initial conditions >> beta = 8/3; >> sigma = 10; >> rho = 28; >> y 0 = 1 e-8 * [ 1 0 0 ]; >> [ tt yy ] = ode 45( @(t, y)lorenz( t, y, beta, rho, sigma), t. Span, y 0 ); >> plot 3( yy(: , 1), yy(: , 2), yy(: , 3), 'b' ) >> title( 'Original' ) Az: -120 El: -20 69

Chaos Now look at sensitivity to initial conditions >> y = yy; >> plot 3( 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… Just pretend it does! Az: -120 El: -20 70

Chaos Now look at sensitivity to initial conditions (1. 01 -1. 00)/1 x 100 >> y 0 = 1 e-8 * [ 1. 01 0 0 ]; >> [ t y ] = ode 45( @(t, y)lorenz( = 1% difference t, y, beta, rho, sigma), t. Span, y 0 ); >> plot 3( yy(: , 1), yy(: , 2), yy(: , 3), 'b', y(: , 1), y(: , 2), y(: , 3), 'y' ) >> title( '1% difference' ) 71

Chaos 1% difference – clearly different paths Az: -120 El: -20 72

Chaos >>y 0=1 e-8*[1. 001 0 0 ]; % 0. 1% difference Az: -120 El: -20 73

Chaos >>y 0=1 e-8*[1. 00001 0 0 ]; % Az: -120 El: -20 0. 001% difference 74

Chaos >>y 0=1 e-8*[1. 0000001 0 0 ]; % 0. 00001% difference Az: -120 El: -20 75

Chaos >>y 0=1 e-8*[1. 0000001 0 0 ]; % 0. 00001% difference y 1(t) vs y 2(t) 76

Chaos >>y 0=1 e-8*[1. 0000001 0 0 ]; % 0. 00001% difference y 1(t) vs y 3(t) 77

Chaos >>y 0=1 e-8*[1. 0000001 0 0 ]; % 0. 00001% difference y 2(t) vs y 3(t) 78

Chaos So even though initial conditions only differ by 1 / 100, 000 of a percent, the trajectories become very different! 79

Chaos 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. 80

Chaos 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 changes 81

Chaos MATLAB is good for studying chaotic systems • Easy to set and change initial conditions or parameters • Solving equations is fast and easy • Plotting and comparing 2 D and 3 D trajectories also fast and easy 82

Chaos Questions? 83

The End 84