Lab08_slides.ipynb

Kernel: SageMath (stable)

Linear System (Revisit)

We have seen system in the following form:

$\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!

$\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 $X_{n+1}$

• Eigenvalues $|\lambda|>1 \implies$ unstable, $|\lambda|\le 1 \implies$ stable

• Longterm behavior converges to dominate eigenvalue's eigenvector

Continuous Time

• Only knows the change $X' = \lim_{\triangle t\rightarrow 0}\frac{X_{n+1}-X_n}{\triangle t}$

• Eigenvalues $real(\lambda)\le 0 \implies \text{ stable, } real(\lambda)>0 \implies \text{ unstable}$

Equilibrium

Most of the time, a continuous time system is written in the following form: $\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?

Jacobian

J = $\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) = x^3+5xy+y^2$

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

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

Recall: Stable/ Unstable Equilibrium Points

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

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

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

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

For $X'=g(X,Y)$, it has slope(gradient) in 2 directions, what should be $<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(\lambda)>0$	$Re(\lambda)=0$	$Re(\lambda)<0$
$Imag(\lambda)=0$	Unstable node	Line of equilibria	stable node
$Imag(\lambda)\ne 0$	Unstable spiral	Center	stable spiral

Why?

What is continuous time system?

$\dot{x} = a x, x(0) = x_0$

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

Continuous time system is written in the following way:

$\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}$

$\begin{bmatrix}x(0)\\y(0)\end{bmatrix} = X_0$

The solution is still:

$X(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) = x_0e^{at}$, if $a<0$, as times goes on, this will decay to 0 (e.g., $a = -1000, t = 100000, x_0e^{-100000000} \rightarrow 0$)

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

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

$\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: $\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: $\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 $\dot{x} =$ linear combination of x and y

But we can try to linearize it:

If we have two points, $(x_1,y_1), (x_2,y_2)$

We know $\dot{x_1} = g(x_1,y_1)$ and $\dot{y_1} = h(x_1,y_1)$

So what is $\dot{x_2}, \dot{y_2}$ ????

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

We can approximate as:

$\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)$

$\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 $(x_1,y_1)$ is equilibrium point!

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

we get:

$\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}$

$\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}$

$\implies$

$\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)$

$\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 = x_2-x_1$, $x_1$ is a constant

$n = y_2-y_1$, $y_1$ is a constant

$\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)$

$\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:

$\dot{\hat{m}} = \frac{\partial g(x_1,y_1)}{\partial x}m+\frac{\partial g(x_1,y_1)}{\partial y}n$

$\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:

$\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(\lambda)>0$	$Re(\lambda)=0$	$Re(\lambda)<0$
$Imag(\lambda)=0$	Unstable node	Line of equilibria	stable node
$Imag(\lambda)\ne 0$	Unstable spiral	Center	stable spiral

Programming

How to get Jacobian??

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

[  2*x + 5*y   5*x + 2*y]
[        6*y 3*y^2 + 6*x]

Tangent Plane

Suppose we have a surface, $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)

fig1+fig2

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

$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 = x^2+3xy+3y^2$ on xz-plane at (2,1,13)

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

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

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

plot(x^2+3*x+3,(x,-4,4))+point((2,13),color='red')

plot(3*y^2+6*y+4,(y,-4,4))+point((1,13),color='red')

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

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

For the tangent plane!

$a(X-X_0)+b(Y-Y_0)= Z^*-Z_0 \implies Z^* = a(X-X_0)+b(Y-Y_0)+Z_0$

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

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

The tangent plane of $z = x^2+3xy+3y^2$ at $(x_0,y_0,z_0)$ is just:

$\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')
fig1+fig2+tangentPlane

Back to Stability

We have $Z = X^3+3XY+X^2Y+3Y^2$

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

$\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)
fig1+fig2

## 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])
show(J)

## Jacobian at equilibrium point
J_now = J.subs({x:0,y:0})
J_now

[0 3]
[3 6]

J_now.eigenvalues()

[-3*sqrt(2) + 3, 3*sqrt(2) + 3]

Unstable Equilibrium Point

plot(x^2,(x,-4,4))

Vector Field

plot_vector_field((3*x^2+2*x*y+3*y,x^2+3*x+6*y),(x,-4,4),(y,-4,4))+point((0,0),color='red',size=30)

Approximate Vector Field with Jabocian

$\begin{bmatrix}\triangle g(X,Y)\\\triangle h(X,Y)\end{bmatrix} = \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}\begin{bmatrix}dx\\dy\end{bmatrix}$

Estimate slope of $(x_2,y_2)$ from $(x_1,y_1)$

$g(X_2,Y_2) = \frac{\partial g(X_1,Y_1)}{\partial X}dx + \frac{\partial g(X_1,Y_1)}{\partial Y}dy+g(X_1,Y_1)$

J_now00 = J_now[0,0]
J_now01 = J_now[0,1]
J_now10 = J_now[1,0]
J_now11 = J_now[1,1]

x_prime(x,y)= 3*x^2+2*x*y+3*y
y_prime(x,y) = x^2+3*x+6*y
x0=0
y0=0

plot_vector_field((J_now00*(x-x0)+J_now01*(y-y0),J_now10*(x-x0)+J_now11*(y-y0)),(x,-1,1),(y,-1,1))+point((0,0),color='red',size=30)