#This script illustrates how to analyze homogeneous linear systems of ODEs in Sage.

#Example System: Consider the homogeneous linear system of ODEs
var('t')
x1 = function('x1')(t);
x2 = function('x2')(t);
de1 = diff(x1,t) == x1 + 3*x2
de2 = diff(x2,t) == 3*x1 + x2

#for functions x1(t) and x2(t), with initial data x1(0) = 2, x2(0) = 6.

#Solution via dsolve: The solution can be obtained with Sage's dsolve_system()
#command:
sol = desolve_system([de1, de2], [x1,x2], ics=[0,2,6]);
print(sol)

#We can pick off the components as
x1sol(t) = x1(t).subs(sol)
x2sol(t) = x2(t).subs(sol)

#Solution via Eigen-analysis: The system above is of the form x' = Ax where
#x = <x1,x2> and
A = matrix(SR,[[1,3],[3,1]]); #Symbolic matrix with real entries
print(A)

#We solve using the eigenvector/value analysis of Section 6.2.2. The
#eigenvalues and vectors of A can be computed as
eigs = A.eigenvectors_right() #Compute right eigenvectors, A*v = lambda*v

#The eigenvalues are
lambda1 = eigs[0][0];
lambda2 = eigs[1][0];
print("The eigenvalues are ", lambda1, " and ", lambda2)

#The eigenvectors are returned in eigs[1][0] and eigs[1][1] and can be
#defined as vectors as
v1 = vector(eigs[0][1][0]);
v2 = vector(eigs[1][1][0]);
print("The eigenvectors are ", v1," and ",v2);
#Sage prints them as row vectors

#Define a matrix P with the eigenvectors as the columns:
P = transpose(matrix([v1, v2]))
print(P)
#as was done in equation (6.18) in the text.

#Let
x0 = vector(SR,[2,6]); #Symbolic vector of initial conditions with real entries

#denote the initial condition vector.

#As per the procedure in Section 6.2.2, we can construct the solution to this ODE
#system by first solving Pc = x0 for c, which in Sage takes the form
c = P\x0;
print(c)

#From (6.17) in the text the solution is then
var('t');
xsol = c[0]*exp(lambda1*t)*v1 + c[1]*exp(lambda2*t)*v2;
print("The solution is x(t) = ",xsol)

#Example 2: Consider the system
var('t')
x1 = function('x1')(t);
x2 = function('x2')(t);
de1 = diff(x1,t) == -5*x1 + 6*x2
de2 = diff(x2,t) == -3*x1 + x2
#for functions x1(t) and x2(t), with initial data x1(0) = 1, x2(0) = 2.

#Solution via Eigen-analysis: The system above is of the form x' = Ax where
#x = <x1,x2> and
A = matrix(SR,[[-5,6],[-3,1]]); #Symbolic matrix with real entries
print(A)

#We solve using the eigenvector/value analysis of Section 6.2.2. The
#eigenvalues and vectors of A can be computed as
eigs = A.eigenvectors_right() #Compute right eigenvectors, A*x = lambda*x

#The eigenvalues are
lambda1 = eigs[0][0];
lambda2 = eigs[1][0];
print("The eigenvalues are ", lambda1, " and ", lambda2)

#The eigenvectors are returned in eigs[1][0] and eigs[1][1] and can be
#defined as vectors as
v1 = vector(eigs[0][1][0]);
v2 = vector(eigs[1][1][0]);
print("The eigenvectors are ", v1," and ",v2);
#Sage prints them as row vectors

#Define a matrix P with the eigenvectors as the columns:
P = transpose(matrix([v1, v2]))
print(P)
#as was done in equation (6.18) in the text.

#Let
x0 = vector(SR,[1,2]); #Symbolic vector with real entries
#denote the initial condition vector.

#As per the procedure in Section 6.2.2, we can construct the solution to this ODE
#system by first solving Pc = x0 for c, which in Sage takes the form
c = P\x0;
print(c)

#From (6.17) in the text the solution is then
var('t');
xsol = c[0]*exp(lambda1*t)*v1 + c[1]*exp(lambda2*t)*v2;
print("The solution is x(t) = ", xsol)

#A real-valued solution can be obtained as follows. First, declare that "t" is a
#real variable with the Sage command
assume(t,'real');

#Now strip off the real part of each component of xsol using Sage's "real" command. This must be done for each
#component, since the "real" command doesn't seem to work on vectors.
xsol1real = simplify(real(xsol[0])); #Strip off the real part of "xsol[0]" (which we know is real)
xsol2real = simplify(real(xsol[1])); #Strip off the real part of "xsol[1]" (which we know is real)
print("The real-valued solution is x1(t) = ", xsol1real);
print("The real-valued solution is x2(t) = ", xsol2real);