352 - Topic 18 - Linear Systems of DEs

8 days ago by Professor352

Solving Linear Differential Equation Systems using Eigenvalues and Eigenvectors

John Travis

Mississippi College

When solving linear systems of differential equations, the collection of equations can be written in the form

$Y^\prime = A Y$

If one assumes that there are solutions that behave like "straight lines" that pass through the origin in the phase plane, then in a vector form the solution must look something like

$Y = r(t) V$

where $r(t)$ is a monotonic function with range $(-\infty,0)$ or $(0,\infty)$ and $V$ is a "direction" vector.

If we presume $r(t) = c e^{\lambda t}$ then

$A (c e^{\lambda t}V) = A Y = Y^\prime = \lambda c e^{\lambda t}V = \lambda Y$

Cancelling $c e^{\lambda t}$ on both sides gives

$AV =  \lambda V $

or

$ 0 = (\lambda I - A)V$

Therefore, the differential equation has a non-trivial straight-line solution of the form $c e^{\lambda t}V$ provided

$\left |{\lambda I - A}\right | = 0$
%hide %auto # Enter the 2x2 matrix corresponding to the desired linear system of DEs... var('t,x,y') A = matrix(QQ,[[-1,2],[0,3]]) #html('For the linear system of differential equations given by') html('$Y^\prime = %s Y$'%str(latex(A))) # Determine eigenvalues and eigenvectors for the given matrix A Eigen = A.eigenmatrix_right() Vecs = Eigen[1].transpose() n=len(Eigen) html('<OL>') for k in range(n): lam = Eigen[0][k][k] V = Vecs[k] soln = e^(lam*t)*V html('<LI> Eigenvalue = $%s$'%str(lam)+' with eigenvector $%s^T$'%str(latex(V))+' yields solution $%s$'%str(latex(soln))) A0 = A[0][0]*x+A[0][1]*y A1 = A[1][0]*x+A[1][1]*y G = plot_vector_field((A0,A1),[x,-2,2],[y,-2,2]) G += line([(0,0),Vecs[0]],thickness=4)+line([(0,0),Vecs[1]],thickness=4) G += line([(0,0),-Vecs[0]],thickness=4)+line([(0,0),-Vecs[1]],thickness=4) G.show(aspect_ratio=true) # plot the straight line solutions using the vectors given. 
       
# Create an interact where the user can adjust components of A and watch the eigensystem (and straight-line solutions) change dynamically 
       
var('t,x,y') # A = matrix(QQ,[[0,1,0,0,0],[0,0,1,0,0],[0,0,0,1,0],[0,0,0,0,1],[31,-19,-12,6,0]]) A = matrix(QQ,[[1,-2],[-3,6]]) Eigen = A.eigenmatrix_right() Vecs = Eigen[1].transpose() n=len(Eigen) for k in range(n): print Eigen[k][k] show(Eigen) show(Vecs) 
        
(7, 0)
(-3, 1/2)


                                
                            

                            
                                

                                    
(7, 0)
(-3, 1/2)


                                

                            
                                

                                    
                                

                        

                    		
                

            

        

        


    




  
    

    
    

      
So, let's solve some linear 2D systems with DISTINCT REAL EIGENVALUES and draw some nice pictures of everything we know.  Then, let's let Sage do the solutions for us and see how things compare.

    





  
    

    
    






    


        
        
        
            
                
            
            
                
var('t,x,y')
# A = matrix(QQ, [[1,1],[4,1]])   # make certain the eigenvalues of this 2x2 matrix are real and distinct
A = matrix(QQ, [[1,2],[4,1]])   # make certain the eigenvalues of this 2x2 matrix are real and distinct

xprime = A[0][0]*x + A[0][1]*y
yprime = A[1][0]*x + A[1][1]*y
show(A)

G = Graphics()

eigs = A.eigenvectors_right()
evals = []
evecs = []
for eig in eigs:
    lamb = eig[0]
    v = eig[1][0]
    show(lamb)
    show(v)
    evals.append(lamb)
    evecs.append(v)
    show(A*v-lamb*v)
    show('----')
    G += plot(3*v,color='yellow')+plot(-3*v,color='cyan')
    G += plot(v)+plot(-v)
    
G += plot_vector_field((xprime,yprime),(x,-4,4),(y,-4,4))

show(G,title='Eigenvectors and vector field')


# Now let's let Sage solve the problem and generate solution curves

X = function('X')(t)
Y = function('Y')(t)

de1 = diff(X,t) - xprime(x=X,y=Y)
de2 = diff(Y,t) - yprime(x=X,y=Y)

t0=0
@interact
def _(x0 = slider(-2,2,1/8,1),y0 = slider(-2,2,1/8,1)):
    H = Graphics()
    sol = desolve_system([de1, de2] , [X,Y], ivar=t, ics=[t0,x0,y0])
    show(sol)
    xcurve = sol[0].rhs()
    ycurve = sol[1].rhs()

    H += parametric_plot((xcurve,ycurve),(t,-1,1),color='red',thickness = 2)
    H += point((x0,y0),color='black',size=20)
    show(H+G, xmin = -5, xmax=5, ymin = -5, ymax=5) 

            
            

            

        

        
        
        

            

                

                    
                               
                    
                    
                        

                            

                                
                                    
                                
                            

                            
                                

                                    
Click to the left again to hide and once more to show the dynamic interactive window

                                

                            
                                

                                    
                                

                        

                    		
                

            

        

        


    




  
    

    
    

      
What about when the eigensystem has imaginary components.  Then, no "straight-line" solutions since $e^{iat} = cos(at)+isin(at)$.  But how does this affect things?

    





  
    

    
    






    


        
        
        
            
                
            
            
                
var('t,x,y')
assume(t,'real')
A = matrix(QQ, [[0,1],[-4,0]])   # make certain the eigenvalues of this 2x2 matrix are complex
# A = matrix(QQ, [[1,1],[-4,0]])   # make certain the eigenvalues of this 2x2 matrix are complex but with a real component
# A = matrix(QQ, [[1,-5],[1,-3]])   # make certain the eigenvalues of this 2x2 matrix are complex but with a real component

xprime = A[0][0]*x + A[0][1]*y
yprime = A[1][0]*x + A[1][1]*y
show(A)

G = Graphics()

eigs = A.eigenvectors_right()
evals = []
evecs = []
solution = 0
for eig in eigs:
    lamb = eig[0].n(digits=3)
    v = eig[1][0].n(digits=3)
    solution += e^(lamb*t)*v
    show(lamb)
    show(v)
    evals.append(lamb)
    evecs.append(v)
    show(A*v-lamb*v)
    show('----')
show(solution)

Y1 = [z.real().simplify() for z in solution]
Y2 = [z.imag().simplify() for z in solution]
show("Stripping off the real part gives:")
show(Y1)
show("Stripping off the imaginary part (w/o i) gives:")
show(Y2)

    
G += plot_vector_field((xprime,yprime),(x,-4,4),(y,-4,4))

show(G,title='Eigenvectors and vector field')


# Now let's let Sage solve the problem and generate solution curves

X = function('X')(t)
Y = function('Y')(t)

de1 = diff(X,t) - xprime(x=X,y=Y)
de2 = diff(Y,t) - yprime(x=X,y=Y)

t0=0
@interact
def _(x0 = slider(-2,2,1/8,1),y0 = slider(-2,2,1/8,1)):
    H = Graphics()
    sol = desolve_system([de1, de2] , [X,Y], ivar=t, ics=[t0,x0,y0])
    show(sol)
    xcurve = sol[0].rhs()
    ycurve = sol[1].rhs()

    H += parametric_plot((xcurve,ycurve),(t,-2*pi,2*pi),color='red',thickness = 2)
    show(H+G, xmin = -5, xmax=5, ymin = -5, ymax=5) 

            
            

            

        

        
        
        

            

                

                    
                               
                    
                    
                        

                            

                                
                                    
                                
                            

                            
                                

                                    
Click to the left again to hide and once more to show the dynamic interactive window

                                

                            
                                

                                    
                                

                        

                    		
                

            

        

        


    




  
    

    
    






    


        
        
        
            
                
            
            
                
A = matrix(QQ, [[0, 1, 0],[0, 0, 1],[5, -7, 3]])   # make certain the eigenvalues of this 2x2 matrix are complex but with a real component

show(A)

eigs = A.eigenvectors_right()
evals = []
evecs = []
for eig in eigs:
    lamb = eig[0]
    v = eig[1][0]
    show(lamb)
    show(v)
    evals.append(lamb)
    evecs.append(v)
    show(A*v-lamb*v)
    show('----') 

            
            

            

        

        
        
        

            

                

                    
                               
                    
                    
                        

                            

                                
                                    

                                
                            

                            
                                

                                    

                                

                            
                                

                                    
                                

                        

                    		
                

            

        

        


    




  
    

    
    






    


        
        
        
            
                
            
            
                
 

            
            

            

        

        
        
        

            

                

                    
                               
                    
                    
                        

                            

                                
                                    
[2, 3]

                                
                            

                            
                                

                                    
[2, 3]