3d_Homogeneous_Constant_Coefficients_Systems.ipynb

Kernel: SageMath 9.2

3d Homogeneous constant coefficient systems

This notebook will assist in creating plots of the vector field and solution curves of 3-dimensional linear homogeneous systems with constant coefficients. These are of the form $$x' = Ax$$ where $A$ is a $3\times 3$ matrix.

Note: we want to input a matrix $A$ as a "list of lists" of the type [[1,2,3],[4,5,6],[7,8,9]] corresponding to the matrix $$\begin{pmatrix} 1&2&3 \\ 4&5&6 \\ 7&8&9\end{pmatrix}.$$

import warnings
warnings.filterwarnings("ignore", category=UserWarning)
@interact
def plot_vector_fields_and_solution_curves(
        a = input_box(default = [[1,1,2],[0,-1,2],[0,0,2]], label = '$A$'),
        num_ics = slider(default = 10, vmin = 0, vmax = 25, label = "# of curves"),
        xrange = slider(1, 10, .5, default = 3), 
        yrange = slider(1,10, .5, default = 3), 
        zrange = slider(1,10, .5, default = 3),
        vector_density = slider(5,30,default = 7, label = "Vector density"),
        auto_update = False):
    x,y,z = var('x,y,z')
    A = Matrix(a).numerical_approx(digits = 2)
    v = A*vector((x,y,z))
    show("A = " + latex(A))
    t = var('t')
    x1,x2,x3 = function('x1')(t), function('x2')(t), function('x3')(t)
    RHS = A*vector((x1,x2,x3))
    LHS = [diff(v,t) for v in [x1,x2,x3]]
    des = [LHS[i] - RHS[i] == 0 for i in [0,1,2]]
    ic_list = [[0,RR.random_element(-xrange/2,xrange/2),
                RR.random_element(-yrange/2,yrange/2),
                RR.random_element(-zrange/2,zrange/2)] for i in range(num_ics)]
    solns = [desolve_system(des, [x1,x2,x3], ics = ic) for ic in ic_list]
    curves = [(sol[0].rhs(), sol[1].rhs(), sol[2].rhs()) for sol in solns]
    step = 0.01
    def evol(curve, positive):
        tcurr = 0
        while -xrange < curve[0](t = tcurr) < xrange and -yrange < curve[1](t = tcurr) < yrange and -zrange < curve[2](t = tcurr) < zrange and abs(tcurr) < 4:
            yield (curve[0](t= tcurr), curve[1](t = tcurr), curve[2](t = tcurr))
            tcurr = tcurr + step if positive else tcurr - step
    lists = [list(reversed(list(evol(c, False)))) + list(evol(c, True)) for c in curves]
    curves_plots = [list_plot(c, plotjoined = True, color = 'red') for c in lists]
    show(sum(curves_plots)  + plot_vector_field3d(v, (x,-xrange,xrange), (y,-yrange,yrange), (z, -zrange, zrange), colors = 'blue', plot_points = vector_density, opacity = 0.4))
    evs = [CDF(l).numerical_approx(digits =2) for l in A.eigenvalues()]
    show(point(evs, size = 12, title = f"Eigenvalues $\\lambda = {evs[0]}, {evs[1]}, {evs[2]}$", figsize = 3, color = 'red'))

# Some matrices to try
[[0,-1,0],[1,0,0],[0,0,1.5]]
[[.5,-1,0],[1,.5,0],[0,0,1.5]]
[[-2,0,0],[0,1,0],[0,0,3]]
[[.5,-1,0],[1,.5,0],[0,0,-1.5]]
[[.5,-1,0],[1,.5,0],[0,0,.5]]
[[2,1,0],[0,2,0],[0,0,1]]

[[0.500000000000000, -1, 0],
 [1, 0.500000000000000, 0],
 [0, 0, -1.50000000000000]]