{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#This script illustrates how to analyze homogeneous linear systems of ODEs in Sage.\n",
    "\n",
    "#Example System: Consider the homogeneous linear system of ODEs\n",
    "var('t')\n",
    "x1 = function('x1')(t);\n",
    "x2 = function('x2')(t);\n",
    "de1 = diff(x1,t) == x1 + 3*x2\n",
    "de2 = diff(x2,t) == 3*x1 + x2\n",
    "#for functions x1(t) and x2(t), with initial data x1(0) = 2, x2(0) = 6."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solution via dsolve: The solution can be obtained with Sage's dsolve_system() command:\n",
    "sol = desolve_system([de1, de2], [x1,x2], ics=[0,2,6]);\n",
    "print(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We can pick off the components as\n",
    "x1sol(t) = x1(t).subs(sol)\n",
    "x2sol(t) = x2(t).subs(sol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solution via Eigen-analysis: The system above is of the form x' = Ax where\n",
    "#x = <x1,x2> and\n",
    "A = matrix(SR,[[1,3],[3,1]]); #Symbolic matrix with real entries\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We solve using the eigenvector/value analysis of Section 6.2.2. The\n",
    "#eigenvalues and vectors of A can be computed as\n",
    "eigs = A.eigenvectors_right() #Compute right eigenvectors, A*v = lambda*v"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The eigenvalues are\n",
    "lambda1 = eigs[0][0];\n",
    "lambda2 = eigs[1][0];\n",
    "print(\"The eigenvalues are \", lambda1, \" and \", lambda2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The eigenvectors are returned in eigs[1][0] and eigs[1][1] and can be\n",
    "#defined as vectors as\n",
    "v1 = vector(eigs[0][1][0]);\n",
    "v2 = vector(eigs[1][1][0]);\n",
    "print(\"The eigenvectors are \", v1,\" and \",v2);\n",
    "#Sage prints them as row vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define a matrix P with the eigenvectors as the columns:\n",
    "P = transpose(matrix([v1, v2]))\n",
    "print(P)\n",
    "#as was done in equation (6.18) in the text."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let\n",
    "x0 = vector(SR,[2,6]); #Symbolic vector of initial conditions with real entries\n",
    "#denote the initial condition vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#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, \n",
    "#which in Sage takes the form\n",
    "c = P\\x0;\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#From (6.17) in the text the solution is then\n",
    "var('t');\n",
    "xsol = c[0]*exp(lambda1*t)*v1 + c[1]*exp(lambda2*t)*v2;\n",
    "print(\"The solution is x(t) = \",xsol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Example 2: Consider the system\n",
    "var('t')\n",
    "x1 = function('x1')(t);\n",
    "x2 = function('x2')(t);\n",
    "de1 = diff(x1,t) == -5*x1 + 6*x2\n",
    "de2 = diff(x2,t) == -3*x1 + x2\n",
    "#for functions x1(t) and x2(t), with initial data x1(0) = 1, x2(0) = 2."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Solution via Eigen-analysis: The system above is of the form x' = Ax where x = <x1,x2> and\n",
    "A = matrix(SR,[[-5,6],[-3,1]]); #Symbolic matrix with real entries\n",
    "print(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#We solve using the eigenvector/value analysis of Section 6.2.2. The eigenvalues and vectors of A can be computed as\n",
    "eigs = A.eigenvectors_right() #Compute right eigenvectors, A*x = lambda*x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The eigenvalues are\n",
    "lambda1 = eigs[0][0];\n",
    "lambda2 = eigs[1][0];\n",
    "print(\"The eigenvalues are \", lambda1, \" and \", lambda2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#The eigenvectors are returned in eigs[1][0] and eigs[1][1] and can be defined as vectors as\n",
    "v1 = vector(eigs[0][1][0]);\n",
    "v2 = vector(eigs[1][1][0]);\n",
    "print(\"The eigenvectors are \", v1,\" and \",v2);\n",
    "#Sage prints them as row vectors"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Define a matrix P with the eigenvectors as the columns:\n",
    "P = transpose(matrix([v1, v2]))\n",
    "print(P)\n",
    "#as was done in equation (6.18) in the text."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Let\n",
    "x0 = vector(SR,[1,2]); #Symbolic vector of initial data with real entries\n",
    "#denote the initial condition vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#As per the procedure in Section 6.2.2, we can construct the solution to this ODE\n",
    "#system by first solving Pc = x0 for c, which in Sage takes the form\n",
    "c = P\\x0;\n",
    "print(c)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#From (6.17) in the text the solution is then\n",
    "var('t');\n",
    "xsol = c[0]*exp(lambda1*t)*v1 + c[1]*exp(lambda2*t)*v2;\n",
    "print(\"The solution is x(t) = \", xsol)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#A real-valued solution can be obtained as follows. First, declare that \"t\" is a real variable with the Sage command\n",
    "assume(t,'real');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Now strip off the real part of each component of xsol using Sage's \"real\" command. This must be done for each\n",
    "#component, since the \"real\" command doesn't seem to work on vectors.\n",
    "xsol1real = simplify(real(xsol[0])); #Strip off the real part of \"xsol[0]\" (which we know is real)\n",
    "xsol2real = simplify(real(xsol[1])); #Strip off the real part of \"xsol[1]\" (which we know is real)\n",
    "print(\"The real-valued solution is x1(t) = \", xsol1real);\n",
    "print(\"The real-valued solution is x2(t) = \", xsol2real);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "SageMath 9.2",
   "language": "sage",
   "name": "sagemath"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}