CoCalc -- Collaborative Calculation in the Cloud
SharedTA Sandbox / Hao's Sandbox / Lab08_slides / Lab08_slides.ipynbOpen in CoCalc

Linear System (Revisit)

We have seen system in the following form:

[Xn+1Yn+1]=M[XnYn]\begin{bmatrix}X_{n+1}\\Y_{n+1}\end{bmatrix} = M \begin{bmatrix}X_n\\Y_n\end{bmatrix}

However, this it sometimes hard to define!!!!

Sometimes we only know the "Propensity" (More tuna, more sharks, too many tuna, not enough food, less tuna)

For these cases, we can use the continuous model!

[XY]=M[XnYn]\begin{bmatrix}X'\\Y'\end{bmatrix} = M \begin{bmatrix}X_n\\Y_n\end{bmatrix}

Continuous Time Model v.s. Discrete Time Model

Discrete Time:

  • Knows the direct value of Xn+1X_{n+1}
  • Eigenvalues λ>1    |\lambda|>1 \implies unstable, λ1    |\lambda|\le 1 \implies stable
  • Longterm behavior converges to dominate eigenvalue's eigenvector

Continuous Time

  • Only knows the change X=limt0Xn+1XntX' = \lim_{\triangle t\rightarrow 0}\frac{X_{n+1}-X_n}{\triangle t}
  • Eigenvalues real(λ)0     stable, real(λ)>0     unstablereal(\lambda)\le 0 \implies \text{ stable, } real(\lambda)>0 \implies \text{ unstable}


Most of the time, a continuous time system is written in the following form: [XY]=[g(X,Y)h(X,Y)]\begin{bmatrix}X'\\Y'\end{bmatrix} = \begin{bmatrix}g(X,Y)\\h(X,Y)\end{bmatrix}

This is not a matrix!!! How can we analyze the stability?


J = [g(X,Y)Xg(X,Y)Yh(X,Y)Xh(X,Y)Y]\begin{bmatrix}\frac{\partial g(X,Y)}{\partial X}&\frac{\partial g(X,Y)}{\partial Y}\\\frac{\partial h(X,Y)}{\partial X}&\frac{\partial h(X,Y)}{\partial Y}\end{bmatrix}

Partial Differential (Only differentiate one of the variable, treat others like constant)

For z(x,y)=x3+5xy+y2z(x,y) = x^3+5xy+y^2

zx=3x2+5y\frac{\partial z}{\partial x} = 3x^2+5y

zy=5x+2y\frac{\partial z}{\partial y} = 5x+2y

Recall: Stable/ Unstable Equilibrium Points

For X=f(X)X' = f(X), which of the following points is stable?

(a) f(x1)=0,f(x1)>0f(x_1)=0, f'(x_1)>0

(b) f(x1)=0,f(x1)<0f(x_1)=0, f'(x_1)<0

(c) f(x1)<0,f(x1)=0f(x_1)<0, f'(x_1)=0

For X=g(X,Y)X'=g(X,Y), it has slope(gradient) in 2 directions, what should be <0<0 ????

Real Parts of Eigenvalues of the Jacobian <0     \implies Stable

Steps of Analyzing Stability

  • Find the equilibrium points
  • Get the Jacobian (symbolic)
  • Calculate the Jacobian at equilibrium points (value)
  • Calculate the Jacobian's eigenvalues!
Re(λ)>0Re(\lambda)>0 Re(λ)=0Re(\lambda)=0 Re(λ)<0Re(\lambda)<0
Imag(λ)=0Imag(\lambda)=0 Unstable node Line of equilibria stable node
Imag(λ)0Imag(\lambda)\ne 0 Unstable spiral Center stable spiral


What is continuous time system?

x˙=ax,x(0)=x0\dot{x} = a x, x(0) = x_0

The solution is x(t)=x0eatx(t) = x_0 e^{at}

Continuous time system is written in the following way:

[x˙y˙]=[a11a12a21a22]A[xy]\begin{bmatrix}\dot{x}\\\dot{y}\end{bmatrix} =\underbrace{\begin{bmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{bmatrix}}_{A}\begin{bmatrix}x\\y\end{bmatrix}

[x(0)y(0)]=X0\begin{bmatrix}x(0)\\y(0)\end{bmatrix} = X_0

The solution is still:

X(t)=eAtX0X(t) = e^{At}X_0 (yes, the power can be a matrix)

What is the criterion that x(t) to be stable?

for x(t)=x0eatx(t) = x_0e^{at}, if a<0a<0, as times goes on, this will decay to 0 (e.g., a=1000,t=100000,x0e1000000000a = -1000, t = 100000, x_0e^{-100000000} \rightarrow 0)

For matrix cases, the term A<0A<0 is about its "Eivenvalues' real parts <0

But this is a linear case!! What if the system is nonlinear?

[x˙y˙]=[g(x,y)h(x,y)]\begin{bmatrix}\dot{x}\\\dot{y}\end{bmatrix} =\begin{bmatrix}g(x,y)\\h(x,y)\end{bmatrix}

Can we still find the eigenvalues of it?

Yes, but no.........

All we want is to write: [x˙y˙]=[a11a12a21a22]A[xy]\begin{bmatrix}\dot{x}\\\dot{y}\end{bmatrix} =\underbrace{\begin{bmatrix}a_{11}&a_{12}\\a_{21}&a_{22}\end{bmatrix}}_{A}\begin{bmatrix}x\\y\end{bmatrix}

From: [x˙y˙]=[g(x,y)h(x,y)]\begin{bmatrix}\dot{x}\\\dot{y}\end{bmatrix} =\begin{bmatrix}g(x,y)\\h(x,y)\end{bmatrix}

In other words, we want to express x˙=\dot{x} = linear combination of x and y

But we can try to linearize it:

If we have two points, (x1,y1),(x2,y2)(x_1,y_1), (x_2,y_2)

We know x1˙=g(x1,y1)\dot{x_1} = g(x_1,y_1) and y1˙=h(x1,y1)\dot{y_1} = h(x_1,y_1)

So what is x2˙,y2˙\dot{x_2}, \dot{y_2} ????

Although we know it is x2˙=g(x2,y2),y2˙=h(x2,y2)\dot{x_2} = g(x_2,y_2), \dot{y_2} = h(x_2,y_2), this is not linear.

We can approximate as:

x^2˙=g(x1,y1)x(x2x1)+g(x1,y1)y(y2y1)+g(x1,y1)\dot{\hat{x}_2} = \frac{\partial g(x_1,y_1)}{\partial x}(x_2-x_1)+\frac{\partial g(x_1,y_1)}{\partial y}(y_2-y_1)+g(x_1,y_1)

y^2˙=h(x1,y1)x(x2x1)+h(x1,y1)y(y2y1)+h(x1,y1)\dot{\hat{y}_2} = \frac{\partial h(x_1,y_1)}{\partial x}(x_2-x_1)+\frac{\partial h(x_1,y_1)}{\partial y}(y_2-y_1)+h(x_1,y_1)

If (x1,y1)(x_1,y_1) is equilibrium point!

At equilibrium point, g(x1,y1)=0,h(x1,y1)=0g(x_1,y_1) =0, h(x_1,y_1)=0

we get:

x^2˙=g(x1,y1)x(x2x1)+g(x1,y1)y(y2y1)+g(x1,y1)=0\dot{\hat{x}_2} = \frac{\partial g(x_1,y_1)}{\partial x}(x_2-x_1)+\frac{\partial g(x_1,y_1)}{\partial y}(y_2-y_1)+\underbrace{g(x_1,y_1)}_{=0}

y^2˙=h(x1,y1)x(x2x1)+h(x1,y1)y(y2y1)+h(x1,y1)=0\dot{\hat{y}_2} = \frac{\partial h(x_1,y_1)}{\partial x}(x_2-x_1)+\frac{\partial h(x_1,y_1)}{\partial y}(y_2-y_1)+\underbrace{h(x_1,y_1)}_{=0}


x^2˙=g(x1,y1)x(x2x1)+g(x1,y1)y(y2y1)\dot{\hat{x}_2} = \frac{\partial g(x_1,y_1)}{\partial x}(x_2-x_1)+\frac{\partial g(x_1,y_1)}{\partial y}(y_2-y_1)

y^2˙=h(x1,y1)x(x2x1)+h(x1,y1)y(y2y1)\dot{\hat{y}_2} = \frac{\partial h(x_1,y_1)}{\partial x}(x_2-x_1)+\frac{\partial h(x_1,y_1)}{\partial y}(y_2-y_1)

Let's define a new coordinate:

m=x2x1m = x_2-x_1, x1x_1 is a constant

n=y2y1n = y_2-y_1, y1y_1 is a constant

x^2˙x^1˙=0=g(x1,y1)x(x2x1)+g(x1,y1)y(y2y1)\dot{\hat{x}_2}-\underbrace{\dot{\hat{x}_1}}_{=0} = \frac{\partial g(x_1,y_1)}{\partial x}(x_2-x_1)+\frac{\partial g(x_1,y_1)}{\partial y}(y_2-y_1)

y^2˙y^1˙=0=h(x1,y1)x(x2x1)+h(x1,y1)y(y2y1)\dot{\hat{y}_2}-\underbrace{\dot{\hat{y}_1}}_{=0} = \frac{\partial h(x_1,y_1)}{\partial x}(x_2-x_1)+\frac{\partial h(x_1,y_1)}{\partial y}(y_2-y_1)

We can now rewrite the above equation as:

m^˙=g(x1,y1)xm+g(x1,y1)yn\dot{\hat{m}} = \frac{\partial g(x_1,y_1)}{\partial x}m+\frac{\partial g(x_1,y_1)}{\partial y}n

n^˙=h(x1,y1)xm+h(x1,y1)yn\dot{\hat{n}} = \frac{\partial h(x_1,y_1)}{\partial x}m+\frac{\partial h(x_1,y_1)}{\partial y}n

or we can use a linear function to represent it:

[m^˙n^˙]=[g(x1,y1)xg(x1,y1)yh(x1,y1)xh(x1,y1)y]jacobian[mn]\begin{bmatrix}\dot{\hat{m}}\\\dot{\hat{n}} \end{bmatrix} = \underbrace{\begin{bmatrix}\frac{\partial g(x_1,y_1)}{\partial x}&\frac{\partial g(x_1,y_1)}{\partial y}\\ \frac{\partial h(x_1,y_1)}{\partial x}&\frac{\partial h(x_1,y_1)}{\partial y}\end{bmatrix}}_{\text{jacobian}}\begin{bmatrix}m\\n\end{bmatrix}

Re(λ)>0Re(\lambda)>0 Re(λ)=0Re(\lambda)=0 Re(λ)<0Re(\lambda)<0
Imag(λ)=0Imag(\lambda)=0 Unstable node Line of equilibria stable node
Imag(λ)0Imag(\lambda)\ne 0 Unstable spiral Center stable spiral


How to get Jacobian??

x = var('x')
y = var('y')
g = x^2+5*x*y+y^2
h = 6*x*y+y^3
[ 2*x + 5*y 5*x + 2*y] [ 6*y 3*y^2 + 6*x]

Tangent Plane

Suppose we have a surface, z=x2+3xy+3y2 z = x^2+3xy+3y^2

x = var('x')
y = var('y')
fig1 = plot3d(x^2+3*x*y+3*y^2,(x,-4,4),(y,-4,4))
fig2 = point3d((2,1,13),color='red',size=20)


Any plane in 3D can be defined by the following function:

a(XX0)+b(YY0)=ZZ0 a(X-X_0)+b(Y-Y_0)= Z-Z_0

How to find the tangent plane that passes through the red point?

1. Define XZ-plane and YZ-plane at that point
2. On both planes, the projection of this tangent plane is a line
3. The slope of this line is the same as the slope of the projection of the surface.

Projection of z=x2+3xy+3y2 z = x^2+3xy+3y^2 on xz-plane at (2,1,13)

=> z=x2+3x×1+3(1)2    z=x2+3x+3z = x^2+3x\times 1+3(1)^2 \implies z=x^2+3x+3

Projection of z=x2+3xy+3y2 z = x^2+3xy+3y^2 on yz-plane at (2,1,13)

=> z=(2)2+3×2y+3y2    z=3y2+6y+4z = (2)^2+3\times 2y+3y^2 \implies z=3y^2+6y+4


Slope on xz plane = d(x2+3x+3)dx=2x+3zx\frac{d(x^2+3x+3)}{dx} = 2x+3 \frac{\partial z}{\partial x}

Slope on yz plane = d(3y2+6y+4)dy=6y+6=zx\frac{d(3y^2+6y+4)}{dy} = 6y+6 = \frac{\partial z}{\partial x}

For the tangent plane!

a(XX0)+b(YY0)=ZZ0    Z=a(XX0)+b(YY0)+Z0 a(X-X_0)+b(Y-Y_0)= Z^*-Z_0 \implies Z^* = a(X-X_0)+b(Y-Y_0)+Z_0

ZX=a=2x+3\frac{\partial Z^*}{\partial X} = a = 2x+3

ZY=b=6y+6\frac{\partial Z^*}{\partial Y} = b = 6y+6

The tangent plane of z=x2+3xy+3y2 z = x^2+3xy+3y^2 at (x0,y0,z0)(x_0,y_0,z_0) is just:

Z(X0,Y0)X(XX0)+Z(X0,Y0)Y(YY0)=ZZ0 \frac{\partial Z(X_0,Y_0)}{\partial X}(X-X_0)+\frac{\partial Z(X_0,Y_0)}{\partial Y}(Y-Y_0)= Z^*-Z_0

x = var('x')
y = var('y')
fig1 = plot3d(x^2+3*x*y+3*y^2,(x,-4,4),(y,-4,4))
fig2 = point3d((2,1,13),color='red',size=20)
tangentPlane = plot3d(7*(x-2)+12*(y-1)+13,(x,-4,4),(y,-4,4),color='green')

Back to Stability

We have Z=X3+3XY+X2Y+3Y2Z = X^3+3XY+X^2Y+3Y^2

ZX=X=3X2+2XY+3Y\frac{\partial Z}{\partial X} = X' = 3X^2+2XY+3Y

ZY=Y=X2+3X+6Y\frac{\partial Z}{\partial Y} = Y' = X^2+3X+6Y

Equilibrium Point: (0,0)

x = var('x')
y = var('y')
fig1 = plot3d(x^3+3*x*y+x^2*y+3*y^2,(x,-4,4),(y,-5,5))
fig2 = point3d((0,0,0),color='red',size=20)
## Get the Jacobian
x = var('x')
y = var('y')
J = jacobian([3*x^2+2*x*y+3*y,x^2+3*x+6*y],[x,y])
(6x+2y2x+32x+36)\left(\begin{array}{rr} 6 \, x + 2 \, y & 2 \, x + 3 \\ 2 \, x + 3 & 6 \end{array}\right)
## Jacobian at equilibrium point
J_now = J.subs({x:0,y:0})
[0 3] [3 6]
[-3*sqrt(2) + 3, 3*sqrt(2) + 3]

Unstable Equilibrium Point