MATH383

MATH383: A first course in differential equationsApril 16, 2020

Homework 12

Due date: April 23, 2020, 11:55PM.

Bibliography: Lessons 22-23

This final homework is a guided tour to some of the behavior encountered in the study of dynamical system.

Problems

Determine the local linearization, eigenvalue maps, and Poincaré sections of the following systems at two different parameter set values of your choice.

1. Duffing oscillator

2. Lorenz system

1Example: Double pendulum

1.1Problem definition

Consider system

u^{'} = v, v^{'} = - \sin u - a v + b \cos (ω t)

that describes the motion of a forced, damped, planar pendulum

\sin u = u - \frac{u^{3}}{3!} + \frac{u^{5}}{5!} - \dots -

1.2Problem analysis

This is a inhomogeneous, nonlinear system of two differential equations of first order,

𝒚^{'} = 𝒇 (𝒚), with 𝒚 = (\begin{array}{l} u \\ v \end{array}), 𝒇 (𝒚) = (\begin{array}{l} v \\ - \sin u - a v + b \cos (ω t) \end{array}) = (\begin{array}{l} f_{1} \\ f_{2} \end{array})

𝑨 = \frac{\partial 𝒇}{\partial 𝒚} = (\begin{array}{ll} \frac{\partial f_{1}}{\partial u} & \frac{\partial f_{1}}{\partial v} \\ \frac{\partial f_{2}}{\partial u} & \frac{\partial f_{2}}{\partial v} \end{array}) = (\begin{array}{ll} 0 & 1 \\ - \cos u & - a \end{array})

p_{𝑨} (λ) = \det (λ 𝑰 - 𝑨) = | \begin{array}{ll} λ & - 1 \\ \cos u & λ + a \end{array} | = λ^{2} + a λ + \cos u = 0 \Rightarrow λ_{1, 2} = \frac{- a \pm \sqrt{a^{2} - 4 \cos^{} u}}{2}

Observations.

If $a = 0$ , then $λ_{1, 2} = \pm i \sqrt{\cos u}$ , locally, around $(u_{k}, v_{k})$ the system behaves as $\cos ((\sqrt{\cos u_{k}}) t)$ , $\sin ((\sqrt{\cos u_{k}}) t)$
If $0 < a < 2$ then the system has complex-conjugate roots and the system spirals towards origin
If $- 2 < a < 0$ then the system has complex-conjugate roots and the system spirals out to infinity
If $a ⩾ 2$ , the system decays towards origin
If $a ⩽ - 2$ , the system evolves towards infinity

1.3System trajectory

1.3.1Trajectory plot

(%i3)

plotdf([v, -sin(u)+0.0*v], [u,v], [trajectory_at,0.7,0.7],[direction,forward],
        [u,-1,1],[v,-1,1],[tstep,0.01],[nsteps,1000])$

(%i4)

1.3.2Numerical compuation of trajectory

(%i81)

(Nsteps: 31, Ncycles: 100, a: 0.2, b: 1.0, w: 2.0)$

(%i82)

[dudt: v, dvdt: -sin(u)-a*v+b*cos(w*t), T: 2*3.14159/w ];

$(%o82) [v, - 0.2 v - \sin (u) + 1.0 \cos (2.0 t), 3.14159]$

(%i83)

[dt: T/Nsteps, tmax: Ncycles*T ];

$(%o83) [0.1013416129032258, 314.159]$

(%i84)

tuvL : rk( [dudt,dvdt], [u,v], [0.8,0.8], [t,0,tmax,dt] )$

(%i85)

N: length(tuvL)$

(%i86)

(%i69)

(%i86)

tuL: makelist( [tuvL[i][1], tuvL[i][2]],i,1,N)$

(%i87)

tvL: makelist( [tuvL[i][1], tuvL[i][3]],i,1,N)$

(%i89)

plot2d( [ [discrete,tuL], [discrete,tvL]], [x,0,tmax], 
        [legend, "u", "v"], [xlabel, "t"])$

(%i90)

1.3.3Numerical computation of phase space plot

(%i90)

uvL: makelist( [tuvL[i][2], tuvL[i][3]],i,1,N)$

(%i91)

plot2d( [ [discrete,uvL] ], [xlabel, "u"], [ylabel, "v"] )$

(%i92)

1.3.4Numerical computation of Poincaré section plot

(%i92)

Nstart: 50*Nsteps$

(%i93)

PoincareL: makelist( uvL[Nsteps*i+Nstart], i, 0, floor((N-Nstart)/Nsteps) )$

(%i94)

plot2d( [discrete,PoincareL], [x,-50,50], [y,-50,50], [style,[points,1,1,1]] )$

(%i95)

(%i31)

2Example: Duffing oscillator

2.1Problem definition

Consider system $u^{''} + δ u^{'} + α u + β u^{3} = γ \sin (ω t)$

u^{'} = v, v^{'} = - δ v - α u - β u^{3} + γ \sin (ω t)

that describes a nonlinear, forced oscillator.

2.2Problem analysis

This is a inhomogeneous, nonlinear system of two differential equations of first order,

𝒚^{'} = 𝒇 (𝒚), with 𝒚 = (\begin{array}{l} u \\ v \end{array}), 𝒇 (𝒚) = (\begin{array}{l} v \\ - δ v - α u - β u^{3} + γ \sin (ω t) \end{array}) = (\begin{array}{l} f_{1} \\ f_{2} \end{array})

𝑨 = \frac{\partial 𝒇}{\partial 𝒚} = (\begin{array}{ll} \frac{\partial f_{1}}{\partial u} & \frac{\partial f_{1}}{\partial v} \\ \frac{\partial f_{2}}{\partial u} & \frac{\partial f_{2}}{\partial v} \end{array}) = (\begin{array}{ll} 0 & 1 \\ - α - 3 β u^{2} & - δ \end{array})

p_{𝑨} (λ) = \det (λ 𝑰 - 𝑨) = | \begin{array}{ll} λ & - 1 \\ α + 3 β u^{2} & λ + δ \end{array} | = λ^{2} + δ λ + (α + 3 β u^{2}) = 0 \Rightarrow λ_{1, 2} = \frac{- δ \pm \sqrt{δ^{2} - 4 (α + 3 β u^{2})}}{2}

Observations.

If $δ = 0$ , then $λ_{1, 2} = \pm \frac{i}{2} \sqrt{α + 3 β u^{2}}$ , locally, around $(u_{k}, v_{k})$ the system behaves as $\cos ((\frac{1}{2} \sqrt{α + 3 β u^{2}}) t)$ , $\sin ((\frac{1}{2} \sqrt{α + 3 β u^{2}}) t)$
If $δ > 0$ the system will decay
If $δ < 0$ the system amplitude will grow exponentially, without bound

2.3System trajectory

2.3.1Trajectory plot

(%i95)

plotdf([v, -0.2*v-0.5*u-0.1*u^3], [u,v], [trajectory_at,0.7,0.7],[direction,forward],
        [u,-1,1],[v,-1,1],[tstep,0.01],[nsteps,1000])$

(%i96)

2.3.2Numerical compuation of trajectory

(%i44)

(Nsteps: 31, Ncycles: 100, a: 1, b: 0.2, d:0.1 , g:1, w: 1.1)$

(%i45)

[dudt: v, dvdt: -d*v-a*u-b*u^3+g*sin(w*t) , T: 2*3.14159/w ];

$(%o45) [v, - 0.1 v - 0.2 u^{3} - u + \sin (1.1 t), 5.711981818181817]$

(%i46)

[dt: T/Nsteps, tmax: Ncycles*T ];

$(%o46) [0.1842574780058651, 571.1981818181818]$

(%i47)

tuvL : rk( [dudt,dvdt], [u,v], [0.1,0.5], [t,0,tmax,dt] )$

(%i48)

N: length(tuvL)$

(%i49)

(%i69)

(%i49)

tuL: makelist( [tuvL[i][1], tuvL[i][2]],i,1,N)$

(%i50)

tvL: makelist( [tuvL[i][1], tuvL[i][3]],i,1,N)$

(%i51)

plot2d( [ [discrete,tuL], [discrete,tvL]], [x,0,tmax], 
        [legend, "u", "v"], [xlabel, "t"])$

(%i52)

2.3.3Numerical computation of phase space plot

(%i52)

uvL: makelist( [tuvL[i][2], tuvL[i][3]],i,1,N)$

(%i53)

plot2d( [ [discrete,uvL] ], [xlabel, "u"], [ylabel, "v"] )$

(%i54)

2.3.4Numerical computation of Poincaré section plot

(%i54)

Nstart: 1*Nsteps$

(%i55)

PoincareL: makelist( uvL[Nsteps*i+Nstart], i, 0, floor((N-Nstart)/Nsteps) )$

(%i58)

plot2d( [discrete,PoincareL], [style,[points,1,1,1]] )$

(%i59)

(%i31)

Figure 1. Poincaré section for Duffing oscillator

3Example: Lorenz system

3.1Problem definition

Consider system

\begin{array}{l} x^{'} = σ (y - x) \\ y^{'} = x (ρ - z) - y \\ z^{'} = x y - β z \end{array}

that describes atmospheric convection (e.g., cumulo-nimbus clouds).

3.2Problem analysis

This is a homogeneous, nonlinear system of three differential equations of first order,

𝒖^{'} = 𝒇 (𝒖), with 𝒖 = (\begin{array}{l} x \\ y \\ z \end{array}), 𝒇 (𝒖) = (\begin{array}{l} σ (y - x) \\ x (ρ - z) - y \\ x y - β z \end{array}) = (\begin{array}{l} f_{1} \\ f_{2} \\ f_{3} \end{array})

𝑨 = \frac{\partial 𝒇}{\partial 𝒖} = (\begin{array}{lll} \frac{\partial f_{1}}{\partial x} & \frac{\partial f_{1}}{\partial y} & \frac{\partial f_{1}}{\partial z} \\ \frac{\partial f_{2}}{\partial x} & \frac{\partial f_{2}}{\partial y} & \frac{\partial f_{2}}{\partial z} \\ \frac{\partial f_{3}}{\partial x} & \frac{\partial f_{3}}{\partial y} & \frac{\partial f_{3}}{\partial z} \end{array}) = (\begin{array}{lll} - σ & σ & 0 \\ ρ - z & - 1 & - x \\ y & x & - β \end{array})

Characteristic polynomial

p_{𝑨} (λ) = \det (λ 𝑰 - 𝑨)

is a cubic, with complicated analytical form of the roots. Analyze behavior around equilibria

At $x = y = z = 0$
$𝑨 = (\begin{array}{lll} - σ & σ & 0 \\ ρ & - 1 & 0 \\ 0 & 0 & - β \end{array}), p_{𝑨} (λ) = \det (λ 𝑰 - 𝑨) = | \begin{array}{lll} λ + σ & - σ & 0 \\ - ρ & λ + 1 & 0 \\ 0 & 0 & λ + β \end{array} | = (λ + β) [λ^{2} + (σ + 1) λ + σ (1 - ρ)]$
has eigenvalues $λ_{1} = - β$ and solutions of quadratic equation $λ^{2} + (σ + 1) λ + σ (1 - ρ) = 0$
$λ_{2, 3} = \frac{- 1 - σ \pm \sqrt{{(σ + 1)}^{2} - 4 σ (1 - ρ)}}{2}$
- If $λ_{1}, λ_{2}, λ_{3} < 0$ the equilibrium point is stable
- If any of $λ_{1}, λ_{2}, λ_{3}$ is positive, equilibrium is unstable.

3.3System trajectory

3.3.1Trajectory plot

(%i60)

plotdf([s*(y-x),x*(r-z)-y], [x,y], [trajectory_at,0.7,0.7],[direction,forward],
       [x,-1,1],[y,-1,1],[tstep,0.01],[nsteps,1000],
       [parameters,"s=0.2,r=0.3,z=0"],
       [sliders,"s=0.1:0.3,r=0.1:0.5,z=-1:1"]
      )$

(%i62)

3.3.2Numerical computation of trajectory

$𝒖^{'} = 𝒇 (𝒖), with 𝒖 = (\begin{array}{l} x \\ y \\ z \end{array}), 𝒇 (𝒖) = (\begin{array}{l} σ (y - x) \\ x (ρ - z) - y \\ x y - β z \end{array}) = (\begin{array}{l} f_{1} \\ f_{2} \\ f_{3} \end{array})$

(%i75)

(Nsteps: 31, Ncycles: 100, s: 0.1, r: 0.2, b:0.1)$

(%i76)

[dxdt: s*(y-x), dydt: x*(r-z)-y, dzdt: x*y-b*z  , T: 3.1415 ];

$(%o76) [0.1 (y - x), x (0.2 - z) - y, x y - 0.1 z, 3.1415]$

(%i77)

[dt: T/Nsteps, tmax: Ncycles*T ];

$(%o77) [0.1013387096774194, 314.15]$

(%i78)

txyzL : rk( [dxdt,dydt,dzdt], [x,y,z], [0.1,0.5,0.2], [t,0,tmax,dt] )$

(%i79)

N: length(txyzL);

$(%o79) 3102$

(%i80)

(%i69)

(%i80)

txL: makelist( [txyzL[i][1], txyzL[i][2]],i,1,N)$

(%i81)

tyL: makelist( [txyzL[i][1], txyzL[i][3]],i,1,N)$

(%i82)

tzL: makelist( [txyzL[i][1], txyzL[i][4]],i,1,N)$

(%i83)

plot2d( [ [discrete,txL], [discrete,tyL], [discrete,tzL] ], 
        [legend, "x", "y", "z"], [xlabel, "t"])$

(%i84)

3.3.3Numerical computation of phase space plot

(%i52)

uvL: makelist( [tuvL[i][2], tuvL[i][3]],i,1,N)$

(%i53)

plot2d( [ [discrete,uvL] ], [xlabel, "u"], [ylabel, "v"] )$

(%i54)

3.3.4Numerical computation of Poincaré section plot

(%i54)

Nstart: 1*Nsteps$

(%i55)

PoincareL: makelist( uvL[Nsteps*i+Nstart], i, 0, floor((N-Nstart)/Nsteps) )$

(%i58)

plot2d( [discrete,PoincareL], [style,[points,1,1,1]] )$

(%i59)

(%i31)

Figure 2. Poincaré section for Duffing oscillator