Report

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., x2 + y2 = r2 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 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π 5 Parametric curves Try It Plot x(t) = cos(t) y(t) = 3sin(t) over 0<t<2π >> ezplot( @(t)cos(t), @(t)3*sin(t) ) x = cos(t), y = 3 sin(t) 2.5 2 1.5 1 y 0.5 0 -0.5 -1 -1.5 -2 -2.5 -3 -2 -1 0 x 1 2 3 6 Parametric curves Try It Plot x(t) = cos3(t) y(t) = 3sin(t) – Hint: use the vector arithmetic operators .*, ./, and .^ to avoid warnings >> ezplot( @(t)cos(t).^3, @(t)3*sin(t) ) x = cos(t)3, y = 3 sin(t) 2.5 2 1.5 1 y 0.5 0 -0.5 -1 -1.5 -2 -2.5 7 -3 -2 -1 0 x 1 2 3 Parametric curves Often parametric curves expressed in polar form ρ = f(θ) Plot with y ρ(θ) θ x ezpolar(fun) ezpolar(f,[thetaMin,thetaMax]) where • f is handle to function ρ = f(θ) • thetaMin, thetaMax specify range of θ – if omitted, range is 0 < θ < 2π 8 Parametric curves Try It Plot ρ = θ >> ezpolar( @(theta)theta ) 90 8 120 60 6 4 150 30 2 180 0 210 330 240 300 r270 = 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),... 90 [0 6*pi ] ) 0.5 120 60 0.4 0.3 150 30 0.2 0.1 180 0 210 330 240 300 270 10 r = sin() cos() Parametric curves 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π 11 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 ) x = cos(4 t), y = sin(4 t), z = t 8 z 6 4 Unfortunately, there's no ezpolar3 2 0 1 0.5 1 0.5 0 0 -0.5 y -0.5 -1 -1 x 12 FOR THRILLS 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 Parametric curves Try It FOR THRILLS Place command window and figure window side by side and use comet3() to plot and CHILLS 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 ) 14 Parametric curves Questions? 15 Phase plane plot For ideal pendulum, θ''+sin( θ(t) ) = 0 Define y1(t) = θ(t), y2(t) = θ'(t) to get y1' (t ) y2 (t ) ' sin y ( t ) y ( t ) 1 2 θ(t) gravity Write pendulum.m function yPrime = pendulum( t, y ) yPrime = [ y(2) -sin( y(1) ) ]'; 16 Phase plane plot 3 different initial conditions θ'(0) θ(0) θ'(0) θ(0) θ(0) θ'(0) y1(0)= θ(0) = 1R=57° y2(0)= θ'(0) = 1R/sec=57°/sec y1(0)= θ(0) = -5R=-286°=74° y2(0)= θ'(0) = 2R/sec=115°/sec y1(0)= θ(0) = 5R=-74° y2(0)= θ'(0) = -2R/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 y1(t) on the horizontal axis and y2(t) on the vertical. Store the results in ta and ya 18 Phase plane plot Qualitatively, what should pendulum do? Try It >> tSpan = [ 0 10 ]; >> y0 = [ 1 1 ]; >> [ ta ya ] = ... ode45( @pendulum, tSpan, y0 ); >> plot( ya(:,1), ya(:,2) ); 1.5 θ(0) θ'(0) y1(0)= θ(0) = 1R=57° y2(0)= θ'(0) = 1R/sec=57°/sec 1 0.5 0 -0.5 -1 19 -1.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 Phase plane plot θ'(0) θ(0) Qualitatively, what should pendulum do? Try It >> tSpan = [ 0 10 ]; >> y0 = [ -5 2 ]; >> [ tb yb ] = ... ode45( @pendulum, tSpan, y0 ); >> plot( yb(:,1), yb(:,2) ); 2.6 y1(0)= θ(0) = -5R=-286°=74° y2(0)= θ'(0) = 2R/sec=115°/sec 2.4 2.2 2 y (t) 2 1.8 1.6 1.4 20 -6 -4 -2 0 2 4 y 1(t) 6 8 10 12 Phase plane plot θ'(0) θ(0) Qualitatively, what should pendulum do? Try It >> tSpan = [ 0 10 ]; >> y0 = [ 5 -2 ]; >> [ tc yc ] = ... ode45( @pendulum, tSpan, y0 ); >> plot( yc(:,1), yc(:,2) ); -1.2 y1(0)= θ(0) = 5R=-74° y2(0)= θ'(0) = -2R/sec=-115°/sec -1.4 -1.6 2 y (t) -1.8 -2 -2.2 -2.4 21 -2.6 -12 -10 -8 -6 -4 -2 y 1(t) 0 2 4 6 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) ] ); 2.5 2 1.5 1 2 y (t) 0.5 0 -0.5 -1 -1.5 -2 -2.5 -5 -4 -3 -2 -1 0 y 1(t) 1 2 3 4 5 22 Phase plane plot 2 1.5 1 0.5 2 y (t) If have initial condition (other than previous 3) that is exactly on curve (red dot) can tell its path in phase plane. 2.5 0 -0.5 -1 -1.5 -2 -2.5 -5 -4 -3 -2 -1 0 y 1(t) 1 2 3 4 5 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) y2(t) y1(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 y1(t) y'(t) = f( t, y ) ! Example – ideal pendulum y1' (t ) y2 (t ) ' sin y ( t ) y ( t ) 1 2 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( 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 8 9 7 8 9 7 8 9 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 y1' (t ) y2 (t ) increments of 0.5 ' sin y (t ) y2 (t ) >> >> >> >> 1 [y1 y2 ] = meshgrid( -5:5, -2.5:0.5:2.5 ); y1Prime = y2; y2Prime = -sin( y1 ); quiver( y1, y2, y1Prime, y2Prime ) 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? 3 2 2 y (t) 1 0 -1 -2 -3 -6 30 -4 -2 0 y 1(t) 2 4 6 Phase plane plot 2.5 2 1.5 1 0.5 2 y (t) 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 0 -0.5 -1 -1.5 -2 -2.5 -5 -4 -3 -2 -1 0 y 1(t) 1 2 3 4 5 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( y1, y2, y1Prime, y2Prime ) >> hold off 32 Phase plane plot 2.5 2 1.5 1 2 y (t) 0.5 0 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. -0.5 -1 -1.5 -2 -2.5 -5 -4 -3 -2 -1 0 y 1(t) 1 2 3 4 33 5 Phase plane plot 2.5 2 1.5 1 2 y (t) 0.5 0 -0.5 -1 -1.5 -2 -2.5 -5 -4 -3 -2 -1 0 y 1(t) 1 2 3 4 5 What path would dot take and why? 34 Phase plane plot 2.5 2 1.5 1 2 y (t) 0.5 0 -0.5 -1 -1.5 -2 -2.5 -5 -4 -3 -2 -1 0 y 1(t) 1 2 3 4 5 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 >> tSpan = [ 0 10 ]; >> y0a = [ 1 1 ]; >> [ ta ya ] = ... ode45( @pendulum, tSpan, y0a ); >> plot( ya(:,1), ya(:,2) ); 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 36 Phase plane plot 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) ) 37 Phase plane plot Initial conditions within 25% 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 2.5 Fairly close 38 Phase plane plot 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) ) 39 Phase plane plot Initial conditions within 10% 2 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 40 Phase plane plot 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) ) 41 Phase plane plot Initial conditions within 1% 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 42 Phase plane plot 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) ) 43 Phase plane plot Initial conditions within 0.1% 1.5 1 0.5 0 -0.5 -1 -1.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 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 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* y (t ) y2 (t ) y1 (t ) ' 1 ' 2 ' 3 y (t ) y1 (t ) y1 (t ) y3 (t ) y2 (t ) y (t ) y1 (t ) y2 (t ) y3 (t ) Notice only two terms have nonlinearities * Also appear in slightly different forms 53 Chaos y (t ) y2 (t ) y1 (t ) ' 1 ' 2 ' 3 y (t ) y1 (t ) y3 (t ) y2 (t ) y (t ) y1 (t ) y2 (t ) y3 (t ) • 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 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); 55 Chaos 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) ) 56 Chaos Try It System converges = 12 20 15 10 5 0 -5 -5 12 Az: -120 El: -20 0 10 8 5 6 4 2 10 0 -2 15 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) ) = 12 = 12 = 12 12 20 20 10 15 15 8 10 3 3 2 y (t) y (t) y (t) 10 6 4 5 5 2 0 0 0 -2 -2 0 2 4 6 y 1(t) 8 10 12 -5 -2 0 2 4 6 y 1(t) 8 10 12 -5 -2 0 2 4 6 8 10 12 y 2(t) 58 Chaos Try It Set ρ = 16 rerun, and do comet3, plot3 30 25 20 15 10 5 0 -20 -5 20 Az: -120 El: -20 -10 15 10 0 5 0 -5 10 -10 -15 20 59 Chaos Try It Notice that "particle" gets pulled over into "hole". "Hole" is called an attractor 30 25 20 15 10 5 0 -20 -5 20 Az: -120 El: -20 -10 15 10 0 5 0 -5 10 -10 -15 20 60 Chaos Try It Set ρ = 20 rerun, and do comet3, plot3 = 20 35 30 25 20 15 10 5 0 -5 20 Az: -120 El: -20 -20 -10 15 0 10 5 0 -5 10 -10 -15 20 61 Chaos Try It Set ρ = 24.2 rerun, and do comet3, plot3 = 24.2 45 40 35 30 25 20 15 10 5 0 -5 25 Az: -120 El: -20 -20 -10 20 15 0 10 5 0 10 -5 -10 -15 20 62 Chaos Try It Set ρ = 24.3 rerun, and do comet3, plot3 Watch comet carefully! = 24.3 45 40 35 30 25 20 15 10 5 0 -5 25 Az: -120 El: -20 -20 -10 20 15 0 10 5 0 -5 10 -10 -15 -20 20 63 Chaos Try It Wow! A small change in ρ causes giant change in trajectory! Particle starts bouncing back and forth between attractors = 24.2 = 24.3 45 45 40 40 35 35 30 30 25 25 20 20 15 15 10 10 5 5 0 -5 25 -20 -10 20 15 0 10 5 0 0 -5 25 -20 -10 20 10 -5 -10 -15 20 15 0 10 5 0 -5 10 -10 -15 -20 20 64 Az: -120 El: -20 Chaos Try It Set ρ = 26 rerun, and do comet3, plot3 = 26 45 40 35 30 25 20 15 10 5 0 -5 30 Az: -120 El: -20 -20 -10 20 0 10 0 -10 10 -20 -30 20 65 Try It Chaos Set ρ = 30 rerun, and do comet3, plot3 In comet, watch hopping back and forth between attractors = 30 60 50 40 30 20 10 0 -20 -10 30 Az: -120 El: -20 -10 0 20 10 10 0 -10 20 -20 -30 30 66 Chaos More common view of Lorenz attractor Has been shown – Bounded (always within a box) – Non-periodic – Never crosses itself = 30 60 50 40 30 20 10 0 -10 -20 -15 -10 -5 0 5 10 15 20 25 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 = 30 60 50 40 30 20 10 0 -10 -20 -15 -10 -5 0 5 10 15 20 25 68 Chaos Now look at sensitivity to initial conditions Original 50 >> 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' ) 40 30 20 10 0 -10 30 -20 -10 20 0 10 0 -10 10 -20 -30 20 Az: -120 El: -20 69 Chaos 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… Just pretend it does! 0% difference 50 40 30 20 10 0 -20 Az: -120 El: -20 -10 30 -10 20 0 10 0 -10 10 -20 -30 20 70 Chaos Now look at sensitivity to initial conditions (1.01-1.00)/1 x 100 >> y0 = 1e-8 * [ 1.01 0 0 ]; >> [ t y ] = ode45( @(t,y)lorenz( = 1% difference 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' ) 71 Chaos 1% difference – clearly different paths 1% difference 50 40 30 20 10 0 -20 -10 30 -10 20 0 10 0 -10 10 -20 -30 Az: -120 El: -20 20 72 Chaos >>y0=1e-8*[1.001 0 0 ]; % 0.1% difference 0.1% difference 50 40 30 20 10 0 -20 -10 30 -10 20 0 10 0 -10 10 -20 -30 Az: -120 El: -20 20 73 Chaos >>y0=1e-8*[1.00001 0 0 ]; % 0.001% difference 0.001% difference 50 40 30 20 10 0 -20 -10 30 -10 20 0 10 0 -10 10 -20 -30 Az: -120 El: -20 20 74 Chaos >>y0=1e-8*[1.0000001 0 0 ]; % 0.00001% difference 0.00001% difference 50 40 30 20 10 0 -20 -10 30 -10 20 0 10 0 -10 10 -20 -30 Az: -120 El: -20 20 75 Chaos >>y0=1e-8*[1.0000001 0 0 ]; % 0.00001% difference y1(t) vs y2(t) 0.00001% difference 30 20 2 y (t) 10 0 -10 -20 -30 -20 -15 -10 -5 0 y 1(t) 5 10 15 20 76 Chaos >>y0=1e-8*[1.0000001 0 0 ]; % 0.00001% difference y1(t) vs y3(t) 0.00001% difference 50 40 3 y (t) 30 20 10 0 -10 -20 -15 -10 -5 0 y 1(t) 5 10 15 20 77 Chaos >>y0=1e-8*[1.0000001 0 0 ]; % 0.00001% difference y2(t) vs y3(t) 0.00001% difference 50 40 3 y (t) 30 20 10 0 -10 -30 -20 -10 0 y 2(t) 10 20 30 78 Chaos So even though initial conditions only differ by 1 / 100,000 of a percent, the trajectories become very different! 0.00001% difference 50 40 3 y (t) 30 20 10 0 -10 -30 -20 -10 0 y 2(t) 10 20 30 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 2D and 3D trajectories also fast and easy 82 Chaos Questions? 83 The End 84