# Reese States That the System of Equations

MATLAB Ordinary Differential Equations – Role Two Greg Reese, Ph. D Inquiry Computing Support Group Academic Engineering Services Miami Academy September 2013

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

Parametric curves Usually describe curve in plane with equation in x and y, e. grand. , 10 2 + y 2 = r ii is circle at origin In other words, tin can write y equally a part of x or vice-versa. Bug • Form may not exist like shooting fish in a barrel to piece of work with • Lots of curves that can’t write as unmarried 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 exist functions of a third variable, say t x = f(t) y = thou(t) Each value of t defines a point (x, y)=( f(t), thou(t) ) that we can plot. Drove of points we get by letting t take on all its values is a parametric curve iv

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 part y(t) • tmin, tmax specify range of t – if omitted, range is 0 < t < 2π 5

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

Parametric curves Effort It Plot x(t) = cos 3(t) y(t) = three 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π eight

Parametric curves Effort 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 half dozen*pi ] ) x

Parametric curves To plot 3 D parametric curves employ 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 role z(t) • tmin, tmax specify range of t – if omitted, range is 0 < t < 2π xi

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

FOR THRILLS Parametric curves Try Information technology 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’ ) xiii

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

Parametric curves Questions? 15

Phase airplane plot For platonic pendulum, θ”+sin( θ(t) ) = 0 Define y one(t) = θ(t), y 2(t) = θ'(t) to get θ(t) gravity Write pendulum. m part y. Prime = pendulum( t, y ) y. Prime = [ y(2) -sin( y(1) ) ]’; 16

Phase aeroplane plot 3 dissimilar initial conditions θ'(0) θ(0) θ'(0) y 1(0)= θ(0) = 1 R=57° y 2(0)= θ'(0) = i R/sec=57°/sec y ane(0)= θ(0) = -5 R=-286°=74° y two(0)= θ'(0) = ii R/sec=115°/sec y i(0)= θ(0) = v R=-74° y 2(0)= θ'(0) = -two R/sec=-115°/sec 17

Phase aeroplane plot Try It For ideal pendulum, θ”+sin( θ(t) ) = 0 solve for the initial conditions θ(0)=1, θ'(0)=1 and fourth dimension = [ 0 10 ] and make a stage plane plot with y ane(t) on the horizontal axis and y ii(t) on the vertical. Store the results in ta and ya 18

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

Phase plane plot θ'(0) θ(0) Qualitatively, what should pendulum exercise? Try Information technology >> t. Bridge = [ 0 10 ]; >> y 0 = [ -five 2 ]; >> [ tb yb ] =. . . ode 45( @pendulum, t. Bridge, y 0 ); >> plot( yb(: , 1), yb(: , 2) ); y one(0)= θ(0) = -5 R=-286°=74° y ii(0)= θ'(0) = 2 R/sec=115°/sec twenty

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

Phase plane plot Try It Graph all iii on one plot >> plot( ya(: , 1), ya(: , two), yb(: , ane), yb(: , 2), . . . yc(: , 1), yc(: , 2) ) >> ax = axis; >> axis( [ -5 5 ax(iii: four) ] ); 22

Phase airplane plot If have initial condition (other than previous three) that is exactly on curve (red dot) tin can tell its path in phase airplane. Q: What if non on curve simply very close to it (xanthous dot)? A: ? 23

Phase airplane plot To help sympathise solution for any initial status, tin can make stage plot and add together information about how each state variable changes with time, i. e. , display the starting time derivative of each state variable. 24

Phase plane plot Volition show charge per unit of alter of state variables at a betoken by cartoon a vector signal there. Horizontal component of vector is rate of modify of variable 1; vertical component of vector is charge per unit of change of variable two. y’ane(t) y’2(t) y 1(t) 25

Phase aeroplane plot Where can we get these rates of change? y (t) y’ane(t) y’two(t) 2 From the state-space conception y one(t) y'(t) = f( t, y ) ! Example – ideal pendulum 26

Phase aeroplane plot To plot vectors at point, use quiver( ten, y, u, five ) u (10, y) v This plots the vectors (u, v) at every point (ten, y) • x is matrix of ten-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 aeroplane plot To brand 10 and y for quiver, use [ x y ] = meshgrid( x. Vec, y. Vec ) Example >> [ ten y ] = meshgrid( 1: 5, 7: 9 ) x = 1 2 3 4 5 y = 7 8 ix seven 8 9 28

Phase plane plot Let’south make quiver plot at every point (x, y) for x going from -v to 5 in increments of ane and y going from -2. v to 2. 5 in increments of 0. five >> >> [y 1 y 2 ] = meshgrid( -five: five, -2. five: 0. 5: two. 5 ); y 1 Prime = y 2; y 2 Prime = -sin( y ane ); quiver( y 1, y two, y ane Prime number, y 2 Prime ) 29

Stage plane plot Now can see charge per unit of change of state variables. MATLAB plots zippo-vectors as a small dot. What is physical meaning of 3 minor dots? 30

Phase plane plot To answer question almost solution with initial atmospheric condition close to those of another solution (yellow dot close to dark-green line), put phase-airplane plot and quiver plot together 31

Phase plane plot Endeavour It >> plot( ya(: , 1), ya(: , ii), yb(: , one), . . . yb(: , two), yc(: , i), yc(: , two) ) >> ax = axis; >> centrality( [ -five 5 ax(iii: four) ] ) >> hold on >> quiver( y one, y 2, y ane Prime, y ii Prime number ) >> hold off 32

Phase aeroplane plot Try It To run across solution path for specific initial conditions, imagine dropping a toy boat (initial status) at a spot in a river (above plot) and watching how current (arrows) pushes it around. 33

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

Phase plane plot From stage-plane plot it appears reasonable to say that if the initial atmospheric condition of the solutions of a differential equation are close to each other, the solutions are likewise shut to each other. 35

Phase plane plot Let’s check out this thought that close initial weather lead to close solutions a petty more. Replot solution to first initial weather >> t. Span = [ 0 10 ]; >> y 0 a = [ i 1 ]; >> [ ta ya ] =. . . ode 45( @pendulum, t. Span, y 0 a ); >> plot( ya(: , 1), ya(: , two) ); 36

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

Phase plane plot Fairly shut 38

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

Phase plane plot twoscore

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

Stage plane plot 42

Phase airplane plot Repeat for 0. 1% greater >> >> n = 0. 001; yy 0 = y 0 a + northward*y 0 a; [ tt yy ] = ode 45( @pendulum, t. Span, yy 0 ); plot( ya(: , i), 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, every bit that famous philosopher might say 46

Phase plane plots Questions? 47

CHAOS or Welcome to my World 48

Chaos theory is co-operative of math that studies behavior of certain kinds of dynamical systems • Chaotic behavior observed in nature, east. grand. , weather • Quantum anarchy theory studies human relationship between chaos and quantum mechanics 49

Anarchy Cluttered systems are: • Deterministic – no randomness involved – If kickoff with identical initial conditions, get identical final states • High sensitivity to initial weather condition – Tiny differences in starting state tin can lead to enormous differences in last land, even over small time ranges • Seemingly random – Unexpected and abrupt changes in land 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 • Besides used equally simple model of laser and dynamo (electric generator) 51

Chaos Gear up of equations • Nonlinear • 3-dimensional • Deterministic, i. east. , no randomness involved Important implications for climate and weather prediction – Atmospheres may exhibit quasi-periodic behavior and may have sharp and seemingly random change, even if fully deterministic – Atmospheric condition can’t exist predicted too far into futurity! 52

Chaos Equations, in state-space class, are* Notice simply ii terms have nonlinearities * Also appear in slightly dissimilar 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 organization Download lorenz. m function y. Prime =. . . lorenz(t, y, beta, rho, sigma) y. Prime = zeros( 3, one ); y. Prime(1) = sigma*( y(two) – y(1) ); y. Prime number(2) = y(1)*( rho – y(3) ) – y(two); y. Prime(iii) = y(one)*y(2) – beta*y(3); 55

Chaos Try It Look at effect of irresolute ρ, starting time with ρ=12 >> beta = eight/3; >> rho = 12; >> sigma = 10; >> t. Span = [ 0 50 ]; >> y 0 = [1 eastward-8 0 0 ]; >> [ t y ] = ode 45( @(t, y)lorenz(. . . t, y, beta, rho, sigma), t. Bridge, y 0 ); >> plot 3( y(: , i), y(: , 2), y(: , iii) ) >> comet 3( y(: , 1), y(: , ii), y(: , 3) ) 56

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

Chaos Endeavour It View two state variables at a time >> plot( y(: , i), y(: , two) ) >> plot( y(: , 1), y(: , 3) ) >> plot( y(: , ii), y(: , 3) ) 58

Anarchy Effort It Set ρ = 16 rerun, and exercise comet three, plot 3 Az: -120 El: -twenty 59

Anarchy Attempt Information technology Notice that “particle” gets pulled over into “hole”. “Hole” is called an attractor Az: -120 El: -20 60

Chaos Try Information technology Set up ρ = 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 three Az: -120 El: -xx 62

Chaos Attempt Information technology Fix ρ = 24. 3 rerun, and do comet 3, plot three Watch comet carefully! Az: -120 El: -20 63

Chaos Attempt It Wow! A small change in ρ causes giant modify in trajectory! Particle starts billowy dorsum and forth between attractors 64 Az: -120 El: -xx

Chaos Try It Set ρ = 26 rerun, and exercise comet 3, plot iii Az: -120 El: -xx 65

Endeavor It Chaos Set ρ = xxx rerun, and do comet three, plot 3 In comet, sentry hopping back and forth between attractors Az: -120 El: -twenty 66

Chaos More common view of Lorenz attractor Has been shown – Bounded (e’er within a box) – Not-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 tin have abrupt changes at seemingly random times • Small variations in a parameter tin produce large changes in trajectories 68

Chaos At present look at sensitivity to initial conditions >> beta = 8/iii; >> sigma = ten; >> rho = 28; >> y 0 = 1 e-viii * [ 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(: , iii), ‘b’ ) >> title( ‘Original’ ) Az: -120 El: -20 69

Chaos Now wait at sensitivity to initial weather >> y = yy; >> plot 3( yy(: , 1), yy(: , ii), yy(: , iii), ‘b’, . . . y(: , 1), y(: , 2), y(: , three), ‘y’ ) >> championship( ‘0% divergence’ ) OOPS! MATLAB issues. Xanthous should exactly cover bluish… Just pretend it does! Az: -120 El: -20 70

Chaos Now look at sensitivity to initial weather (1. 01 -one. 00)/1 x 100 >> y 0 = 1 e-eight * [ one. 01 0 0 ]; >> [ t y ] = ode 45( @(t, y)lorenz( = i% difference t, y, beta, rho, sigma), t. Span, y 0 ); >> plot 3( yy(: , ane), yy(: , 2), yy(: , 3), ‘b’, y(: , one), y(: , two), y(: , 3), ‘y’ ) >> title( ‘1% divergence’ ) 71

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

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

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

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

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

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

Chaos >>y 0=one east-8*[i. 0000001 0 0 ]; % 0. 00001% difference y 2(t) vs y iii(t) 78

Anarchy And then fifty-fifty though initial conditions only differ by one / 100, 000 of a percent, the trajectories become very different! 79

Chaos This farthermost sensitivity to initial atmospheric condition is often called The Butterfly Consequence A butterfly flapping its wings in Brazil can cause a tornado in Texas. eighty

Chaos Lorenz equations are good case of chaotic system • Deterministic • Loftier sensitivity to initial conditions – Very tiny differences in starting state can lead to substantial differences in final country • Unexpected and precipitous changes in country • Sensitive to slight parameter changes 81

Chaos MATLAB is proficient for studying chaotic systems • Easy to gear up and change initial weather condition or parameters • Solving equations is fast and piece of cake • Plotting and comparing 2 D and 3 D trajectories also fast and easy 82

Anarchy Questions? 83

The End 84

## Reese States That the System of Equations

Source: https://slidetodoc.com/matlab-ordinary-differential-equations-part-ii-greg-reese/