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 II © 2010 -2013 Greg Reese. All rights

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.

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

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,

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 Try Information technology Plot 10(t) = cos(t) y(t) = 3 sin(t) over 0<t<2π

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 Try It Plot x(t) = cos 3(t) y(t) = 3 sin(t) –

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

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 Try It Plot ρ = θ >> ezpolar( @(theta)theta ) 9 “>
        </p>
<p class=Parametric curves Effort It Plot ρ = θ >> ezpolar( @(theta)theta ) 9

Parametric curves Try It Plot ρ = sin(θ)cos(θ) over 0 < θ < 6π

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 use ezplot(funx, funy, funz) ezplot(funx, funy,

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 x(t) = cos(4 t) y(t) = sin(4 t) z(t)

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 It Plot x(t) = cos(4 t) y(t) = sin(4

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 command window and figure window

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

Parametric curves Questions? 15

Phase plane plot For ideal pendulum, θ''+sin( θ(t) ) = 0 Define y 1(t)

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 plane plot 3 different initial conditions θ'(0) θ(0) θ'(0) y 1(0)= θ(0) =

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 plane plot Try It For ideal pendulum, θ''+sin( θ(t) ) = 0 solve

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

Phase plane plot θ'(0) θ(0) Qualitatively, what should pendulum do? Try It >> t.”>
        </p>
<p class=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.”>
        </p>
<p class=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 three on one plot >> plot( ya(:”>
        </p>
<p class=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 plane plot If have initial condition (other than previous 3) that is exactly

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 plane plot To help understand solution for any initial condition, can make phase

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 Will show rate of change of state variables at a point

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 plane plot Where can we get these rates of change? y (t) y'1(t)

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

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 plane plot To make x and y for quiver, use [ x y

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's make quiver plot at every point (x, y) for x

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

Phase plane plot Now can see rate of change of state variables. MATLAB plots

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 about solution with initial conditions close to those

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 Try It >> plot( ya(: , 1), ya(: , 2), yb(:”>
        </p>
<p class=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 plane plot Try It To see solution path for specific initial conditions, imagine

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 take and why? 34

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

Phase plane plot From phase-plane plot it appears reasonable to say that if the

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 Now let's solve again with initial conditions 25% greater and plot

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 close 38

Phase plane plot Fairly shut 38

Phase plane plot Repeat for 10% greater >> >> n = 0. 1; yy”>
        </p>
<p class=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 40

Phase plane plot twoscore

Phase plane plot Repeat for 1% greater >> >> n = 0. 01; yy”>
        </p>
<p class=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

Phase plane plot 42

Stage plane plot 42

Phase plane plot Repeat for 0. 1% greater >> >> n = 0. 001;”>
        </p>
<p class=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 44

Phase plane plot Again, it appears that if the initial conditions of the solutions

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 plot Well, every bit that famous philosopher might say 46

Phase plane plots Questions? 47

Phase plane plots Questions? 47

CHAOS or Welcome to my World 48

CHAOS or Welcome to my World 48

Chaos theory is branch of math that studies behavior of certain kinds of dynamical

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

Chaos Chaotic systems are: • Deterministic – no randomness involved – If start with

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

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 Set of equations • Nonlinear • Three-dimensional • Deterministic, i. e. , no

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 form, are* Notice only two terms have nonlinearities * Also

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”>
        </p>
<p class=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 =. .

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 changing ρ, start with ρ=12 >> beta”>
        </p>
<p class=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 Try It System converges Az: -120 El: -20 57

Chaos Try It View two state variables at a time >> plot( y(: ,”>
        </p>
<p class=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

Chaos Try It Set ρ = 16 rerun, and do comet 3, plot 3

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

Chaos Try It Set ρ = 20 rerun, and do comet 3, plot 3

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

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

Chaos Try It Set ρ = 24. 3 rerun, and do comet 3, plot

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

Chaos Try It Wow! A small change in ρ causes giant change in trajectory!

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 do comet 3, plot 3

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

Try It Chaos Set ρ = 30 rerun, and do comet 3, plot 3

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 (always within

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

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 Now look at sensitivity to initial conditions >> beta = 8/3; >> sigma”>
        </p>
<p class=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 look at sensitivity to initial conditions >> y = yy; >> plot”>
        </p>
<p class=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 conditions (1. 01 -1. 00)/1 x 100

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 – clearly different paths Az: -120 El: -20 72

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

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

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

Chaos >>y 0=1 e-8*[1. 0000001 0 0 ]; % 0. 00001% difference Az: -120″>
        </p>
<p class=Chaos >>y 0=1 east-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)”>
        </p>
<p class=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 e-8*[1. 0000001 0 0 ]; % 0. 00001% difference y 1(t)”>
        </p>
<p class=Chaos >>y 0=1 east-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)”>
        </p>
<p class=Chaos >>y 0=one east-8*[i. 0000001 0 0 ]; % 0. 00001% difference y 2(t) vs y iii(t) 78

Chaos So even though initial conditions only differ by 1 / 100, 000 of

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 extreme sensitivity to initial conditions is often called The Butterfly Effect A

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 example of chaotic system • Deterministic • High sensitivity

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 good for studying chaotic systems • Easy to set and change

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

Chaos Questions? 83

Anarchy Questions? 83

The End 84

The End 84

Reese States That the System of Equations

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