StateSpaceModel.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Function Collection to help solve linear systems problems"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import sympy as sy\n",
    "from sympy.integrals import inverse_laplace_transform as ilt\n",
    "from sympy.integrals import laplace_transform as lt\n",
    "from sympy import pretty_print as pp\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Matrix Functions\n",
    "\n",
    "def round_expr(expr, num_digits):\n",
    "    # https://www.thetopsites.net/article/54462201.shtml\n",
    "    return expr.xreplace({n : round(n, num_digits) for n in expr.atoms(sy.Number)})\n",
    "\n",
    "def matrix_info(matrix):\n",
    "    ''' displays the determenant, inverse and charactersitic funtion of a matrix'''\n",
    "    sMatrix = sy.Matrix(matrix)\n",
    "    # char func is det(A- (Lambda x I)) = 0\n",
    "    L = sy.symbols('L')\n",
    "    LI = L * sy.eye(sMatrix.shape[0])\n",
    "    charMatrix = sMatrix - LI\n",
    "    print(\"A-Lamda x I: \")\n",
    "    display(charMatrix)\n",
    "    inverse = None\n",
    "    if sy.det(sMatrix)!=0:\n",
    "        inverse = sMatrix.inv()\n",
    "    else:\n",
    "        inverse = \"Not invertable, det=0\"\n",
    "    display({\n",
    "        \"det\": sMatrix.det(),\n",
    "        \"inverse\":inverse,\n",
    "        \"charFunction\": sy.latex(charMatrix.det())\n",
    "    })\n",
    "    \n",
    "def get_eign(matrix):\n",
    "    '''displays eignvalues and vectors of a matrix'''\n",
    "    sMatrix = sy.Matrix(matrix)\n",
    "    sym_eignvals = list(sMatrix.eigenvals().keys())\n",
    "    sym_eignvects = [list(tup[2][0]) for tup in sMatrix.eigenvects()]\n",
    "    display(\"Sympy eigenvalues: \", sym_eignvals)\n",
    "    display(\"Sympy eigenvectors: \", sym_eignvects)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# State space model functions\n",
    "\n",
    "    \n",
    "def state_trans_matrix_inverse_laplace_method(a):\n",
    "    '''warning can be slow, displays the state transisiton matrix assuming it exisits'''\n",
    "    a = sy.Matrix(a)\n",
    "    assert a.shape[0] == a.shape[1], \"Must be a square matrix\"\n",
    "    s, t = sy.symbols('s t')\n",
    "    invert = s*sy.eye(a.shape[0])-a\n",
    "    display(\"We will now invert this\", invert)\n",
    "    invert = invert.inv()\n",
    "    display(\"After we inverted\", invert)\n",
    "    out = ilt(invert, s, t)\n",
    "#     out = sy.simplify(invert)\n",
    "    display(\"Now we apply inverse Laplace transform to get the state transition matrix\", out)\n",
    "\n",
    "def state_trans_cayley_hamilton_method(a):\n",
    "    '''Only implemented for a 2 by 2'''\n",
    "    a = sy.Matrix(a)\n",
    "    assert a.shape[0]==2 , \"implemented for 2by2 only\"\n",
    "    eignvals = list(a.eigenvals().keys())\n",
    "    t, a0, a1 = sy.symbols('t a0 a1')\n",
    "    eq1 = sy.Eq(sy.exp(t*eignvals[0]), a0+a1*eignvals[0])\n",
    "    eq2 = sy.Eq(sy.exp(t*eignvals[1]), a0+a1*eignvals[1])\n",
    "    display(eq1,eq2)\n",
    "    # Now we make the system of two variables into a solvable thing\n",
    "    # https://stackoverflow.com/questions/31547657/how-can-i-solve-system-of-linear-equations-in-sympy\n",
    "\n",
    "    e1 = sy.exp(t*eignvals[0]) -1*(a0+a1*eignvals[0])\n",
    "    e2 = sy.exp(t*eignvals[1]) -1*(a0+a1*eignvals[1])\n",
    "    sol = sy.solve([e1,e2], (a0,a1))\n",
    "    display(\"a0 and a1 aka alpha 0 and alpha 1\")\n",
    "    display(sol[a0],sol[a1])\n",
    "    ans = sol[a0] * sy.eye(2) + sol[a1]*a\n",
    "    sy.simplify(ans)\n",
    "    display(ans)\n",
    "    \n",
    "\n",
    "def state_trans_matrix(a,symbol = sy.symbols('t'), show_steps = True):\n",
    "    a = sy.Matrix(a)\n",
    "    state_trans = sy.exp(a*symbol)\n",
    "    if show_steps:\n",
    "        display(\"the state trans matrix is\", state_trans)\n",
    "    return state_trans\n",
    "\n",
    "def state_var_homo_response(A_matrix, initial_condition_col,discrete = False, show_steps = True):\n",
    "    '''x_homo at time t = state_trans_matrix * initial_cond '''\n",
    "    A_matrix = sy.Matrix(A_matrix)\n",
    "    initial_condition_col = sy.Matrix(initial_condition_col)\n",
    "    if discrete:\n",
    "        state_trans_mat = state_trans_matrix_discrete(A_matrix, show_steps = show_steps)\n",
    "    else:\n",
    "        state_trans_mat = state_trans_matrix(A_matrix, show_steps=show_steps)\n",
    "        \n",
    "    x_homo =  state_trans_mat * initial_condition_col\n",
    "    if show_steps:\n",
    "        display(\"Finding the homogenous response\")\n",
    "        display(\"A matrix\",a,\"Initial condition\",initial_condition_col )\n",
    "        display(\"x homogenous response at time t is \", x_homo)\n",
    "    return x_homo\n",
    "\n",
    "def state_var_forced_response(a, b, u, show_steps = True):\n",
    "    '''Takes the A matrix and the B inputs column and returns the forced \n",
    "    response of the state space.\n",
    "    '''\n",
    "    a = sy.Matrix(a)\n",
    "    b = sy.Matrix(b)\n",
    "    # T is for tau\n",
    "    t, T = sy.symbols('t T')\n",
    "    state_trans_neg = state_trans_matrix(-1*a, T, show_steps = False)\n",
    "    state_trans = state_trans_matrix(a, t, show_steps = False)\n",
    "    # This is in case the u param is dependant on t, if it is we need to make it\n",
    "    # into a T for the integration below\n",
    "    try:\n",
    "        u = u.replace(t,T)\n",
    "    except:\n",
    "        pass\n",
    "    unevaluated_conv_integral = sy.Integral(state_trans_neg * b*u ,(T,0,t))\n",
    "    conv_integral = unevaluated_conv_integral.doit() \n",
    "    forced_resp = state_trans * conv_integral\n",
    "    if show_steps:\n",
    "        display(\"Finding the forced response\")\n",
    "        display(\"A matrix\",a, \"B matrix\",b)\n",
    "        display(\"State transition matrix: \",state_trans)\n",
    "        display(\"Negative state trans matrix\", state_trans_neg)\n",
    "        display(\"The convolution integral\", unevaluated_conv_integral)\n",
    "        display(\"Evaluated convolution integral\", conv_integral)\n",
    "        display(\"Foreced response is \", forced_resp)\n",
    "        \n",
    "    return forced_resp\n",
    "\n",
    "def output_response(a,b,c,d,u, initial_condition_col, discrete = False, show_steps=True):\n",
    "    '''a,b,c,d,u,ic matrices, and a bool if discrete'''\n",
    "    a = sy.Matrix(a)\n",
    "    b = sy.Matrix(b)\n",
    "    c = sy.Matrix(c)\n",
    "    d = sy.Matrix(d)\n",
    "    initial_condition_col = sy.Matrix(initial_condition_col)\n",
    "    # set show steps to false because it will be redundant\n",
    "    homo_resp = state_var_homo_response(a,initial_condition_col,discrete = discrete, show_steps = False)\n",
    "    forced_resp = []\n",
    "    \n",
    "    if discrete:\n",
    "        forced_resp = state_var_forced_response_discrete(a, b,u, show_steps = False)\n",
    "    else:\n",
    "        forced_resp = state_var_forced_response(a, b,u, show_steps = False)\n",
    "\n",
    "    full_resp_times_c_col = c * forced_resp + c * homo_resp\n",
    "    \n",
    "    y =  full_resp_times_c_col + d*u\n",
    "    if show_steps:\n",
    "        display(\"finding the output response\")\n",
    "        display(\"A matrix\",a, \"B matrix\",b, \"C matrix\",c,\"D matrix\",d)\n",
    "        display(\"initial condition\", initial_condition_col)\n",
    "        display(\"homogenous resp\", homo_resp)\n",
    "        display(\"forced resp\", forced_resp)\n",
    "        display(\"Full x response ie x homogenous + x forced\", homo_resp + forced_resp)\n",
    "        display(\"Full response times c col at time t (or k if discrete)\", full_resp_times_c_col)\n",
    "        display(\"the output y is \", y)\n",
    "    \n",
    "    return y\n",
    "\n",
    "# def neat_output_resp_discrete():\n",
    "#     sy.collect(round_expr(output_response(a,b,c,d,u,ic,discrete=True, show_steps=False),2)[0], k)\n",
    "    "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'We will now invert this'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}s - 1 & -1 & 0\\\\0 & s - 2 & -1\\\\0 & 0 & s - 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[s - 1,    -1,     0],\n",
       "[    0, s - 2,    -1],\n",
       "[    0,     0, s - 1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'After we inverted'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{1}{s - 1} & \\frac{1}{s^{2} - 3 s + 2} & \\frac{1}{s^{3} - 4 s^{2} + 5 s - 2}\\\\0 & \\frac{1}{s - 2} & \\frac{1}{s^{2} - 3 s + 2}\\\\0 & 0 & \\frac{1}{s - 1}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1/(s - 1), 1/(s**2 - 3*s + 2), 1/(s**3 - 4*s**2 + 5*s - 2)],\n",
       "[        0,          1/(s - 2),          1/(s**2 - 3*s + 2)],\n",
       "[        0,                  0,                   1/(s - 1)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Now we apply inverse Laplace transform to get the state transition matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}e^{t} \\theta\\left(t\\right) & \\left(e^{t} - 1\\right) e^{t} \\theta\\left(t\\right) & \\left(- t + e^{t} - 1\\right) e^{t} \\theta\\left(t\\right)\\\\0 & e^{2 t} \\theta\\left(t\\right) & \\left(e^{t} - 1\\right) e^{t} \\theta\\left(t\\right)\\\\0 & 0 & e^{t} \\theta\\left(t\\right)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[exp(t)*Heaviside(t), (exp(t) - 1)*exp(t)*Heaviside(t), (-t + exp(t) - 1)*exp(t)*Heaviside(t)],\n",
       "[                  0,            exp(2*t)*Heaviside(t),      (exp(t) - 1)*exp(t)*Heaviside(t)],\n",
       "[                  0,                                0,                   exp(t)*Heaviside(t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "None"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "a = sy.Matrix([\n",
    "    [1,1,0],\n",
    "    [0,2,1],\n",
    "    [0,0,1]\n",
    "])\n",
    "display(state_trans_matrix_inverse_laplace_method(a))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def trans_matrix(a,b,c,d):\n",
    "    '''returns the transfer matrix, which is the generlization of the transfer function'''\n",
    "    a = sy.Matrix(a)\n",
    "    b = sy.Matrix(b)\n",
    "    c = sy.Matrix(c)\n",
    "    d = sy.Matrix(d)\n",
    "    assert a.shape[0] == a.shape[1], \"A Must be a square matrix\"\n",
    "    s, t = sy.symbols('s t')\n",
    "    invert = s*sy.eye(a.shape[0])-a\n",
    "    display(\"We will now invert this\", invert)\n",
    "    invert = invert.inv()\n",
    "    display(\"After we inverted\", invert)\n",
    "    G = c * invert * b \n",
    "    display(\"Answer of c * invert * b \", G)\n",
    "    G = G +d\n",
    "    display(\"Answer with +d so c * invert * b + d\", G)\n",
    "    G = sy.simplify(sy.together(G))\n",
    "    display(\"simplified answer\", G)\n",
    "    return G\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'the state trans matrix is'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 e^{- t} - e^{- 2 t} & e^{- t} - e^{- 2 t}\\\\- 2 e^{- t} + 2 e^{- 2 t} & - e^{- t} + 2 e^{- 2 t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[   2*exp(-t) - exp(-2*t),    exp(-t) - exp(-2*t)],\n",
       "[-2*exp(-t) + 2*exp(-2*t), -exp(-t) + 2*exp(-2*t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 e^{- t} - e^{- 2 t} & e^{- t} - e^{- 2 t}\\\\- 2 e^{- t} + 2 e^{- 2 t} & - e^{- t} + 2 e^{- 2 t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[   2*exp(-t) - exp(-2*t),    exp(-t) - exp(-2*t)],\n",
       "[-2*exp(-t) + 2*exp(-2*t), -exp(-t) + 2*exp(-2*t)]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = sy.Matrix([\n",
    "    [0,1],\n",
    "    [-2,-3]\n",
    "])\n",
    "state_trans_matrix(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{3 \\delta\\left(0\\right)}{z^{2}}$"
      ],
      "text/plain": [
       "3*DiracDelta(0)/z**2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "with the other functino\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "[-3 + 2*sqrt(6), -2*sqrt(6) - 3]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[5, 3 - 2*sqrt(6), 3 + 2*sqrt(6)]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "def z_transform(expr,n_symbol, start=0, stop=100):\n",
    "    '''Uses the z transform defination to get its summation, use \".doit()\" on the return value\n",
    "    to evaluate the summation\n",
    "    '''\n",
    "    z,k = sy.symbols('z k')\n",
    "    return sy.Sum(expr.replace(n_symbol, k)*z**(-k),(k,start,stop))\n",
    "n,z = sy.symbols('n,z')\n",
    "display(z_transform(3*sy.DiracDelta(n-2), n ,0, 100).doit())\n",
    "a = sy.Matrix([\n",
    "    [1,2,4],[2,5,0],[4,0,5]\n",
    "])\n",
    "# state_trans_matrix_inverse_laplace_method(a)\n",
    "print(\"with the other functino\")\n",
    "display( sy.solve(z**2+6*z-15,z))\n",
    "display( sy.solve(z**3-11*z**2+15*z+75,z))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "n,z = sy.symbols(\"n z\")\n",
    "\n",
    "# n = 1/4 * z**(-1)\n",
    "# d = (1-0.5*z**(-1))*(1-0.25*z**(-1))\n",
    "# n = z**2+5*z+6\n",
    "# d = z+2\n",
    "\n",
    "# n= z\n",
    "# d = z-0.5\n",
    "# sy.div(n,d)\n",
    "# expr = z/((z+1)*(z+2))\n",
    "# expr = n/d\n",
    "# needs edge cases check, what happens when expr is just a number?\n",
    "def inverse_Z_transform(expr, z_symbol = sy.symbols('z'), n_symbol= sy.symbols('n')):\n",
    "    '''Uses the residue theorm explained here http://course.ece.cmu.edu/~ece491/lectures/L07/IZT_ContourIntegration.pdf'''\n",
    "#     import ipdb; ipdb.set_trace()\n",
    "    expr = sy.simplify(expr)\n",
    "    f = sy.simplify(expr * z_symbol**(n_symbol-1))\n",
    "    numer,denom = f.as_numer_denom()\n",
    "    poles = sy.roots(denom)\n",
    "    inverse = 0\n",
    "    display(f)\n",
    "    \n",
    "    # looping through poles to find residues, I should use limits, but I am replacing instead\n",
    "    assert len(poles) > 0, \"There must be at least one pole, or answer is just a delta dirac\"\n",
    "    for p, multiplicty in poles.items():\n",
    "        derivative_order = multiplicty - 1\n",
    "        s = multiplicty\n",
    "        res = (z_symbol-p)**s * f\n",
    "        display(res)\n",
    "        for i in range(derivative_order):\n",
    "            res = sy.Derivative(res, z_symbol)\n",
    "        display(res)\n",
    "        factorial_div = sy.Rational(1,sy.factorial(s-1))\n",
    "        res = sy.Mul(factorial_div, res)\n",
    "        display(\"The residue equation looks like this\", res)\n",
    "        res = res.doit()\n",
    "        res = res.replace(z_symbol,p)\n",
    "        inverse += sy.simplify(res)\n",
    "        display(f\"The residue found for pole {p} is: \",sy.simplify(res))\n",
    "    return inverse\n",
    "\n",
    "def inverse_Z_transform_partial_fraction_method(expr, z_symbol = sy.symbols('z'), n_symbol= sy.symbols('n'), show_steps = True):\n",
    "    expr_over_z = expr/z_symbol\n",
    "    fracs = expr_over_z.apart()\n",
    "    to_transform = sy.apart(sy.simplify(fracs * z_symbol))\n",
    "    if show_steps:\n",
    "        X = sy.Function(\"X\")\n",
    "        display(\"We will take the partial fraction of this\",sy.Eq(X(z_symbol)/z_symbol, expr_over_z))\n",
    "        display(sy.Eq(X(z_symbol)/z_symbol, fracs))\n",
    "        display(\"Now we will take the Z-inverse of this, note that thepoles might have changes\",sy.Eq(X(z_symbol), to_transform))\n",
    "    return to_transform\n",
    "\n",
    "def state_trans_matrix_Z_transform_method(a, show_steps = True):\n",
    "    '''For discrete time'''\n",
    "    a = sy.Matrix(a)\n",
    "    assert a.shape[0] == a.shape[1], \"Must be a square matrix\"\n",
    "    z, k = sy.symbols('z k')\n",
    "    invert_me = z*sy.eye(a.shape[0])-a\n",
    "    inverted = invert_me.inv()\n",
    "#     out = ilt(inverted, z, k)\n",
    "#      out = inverted.applyfunc(inverse_Z_transform)\n",
    "    if show_steps:\n",
    "        display(\"We will now invert this\", invert_me)\n",
    "        print(invert_me)\n",
    "        display(\"After we inverted\", inverted)\n",
    "#         out = inverted.applyfunc(inverse_Z_transform_partial_fraction_method)\n",
    "\n",
    "#         display(\"Now we apply inverse Laplace transform to get the state transition matrix\", out)\n",
    " \n",
    "    return out\n",
    "\n",
    "# v = (z**2+z)/((z-1)**2)\n",
    "# v = z/((z+1)*(z+2))\n",
    "# v = (z**2+2*z)/(z**2-3*z+2) # got the wrong answer for this one\n",
    "# inverse_Z_transform(v)\n",
    "a= sy.Matrix([\n",
    "[sy.Rational(-3,4) , sy.Rational(-1,4)],\n",
    "    [sy.Rational(-1,4) , sy.Rational(-3,4)]\n",
    "])\n",
    "b = [\n",
    "    [1],\n",
    "    [0]\n",
    "]\n",
    "u=1\n",
    "c = [[1,0]]\n",
    "d = [4]\n",
    "ic = [[0],\n",
    "     [0]]\n",
    "\n",
    "\n",
    "\n",
    "a = sy.Matrix([\n",
    "    [0,1],\n",
    "    [sy.Rational(-1,4), -1]\n",
    "])\n",
    "\n",
    "# state_trans_matrix_Z_transform_method(a)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{4 z + 4}{4 z^{2} + 4 z + 1} & - \\frac{1}{4 \\left(- \\frac{z^{2}}{4} - \\frac{z}{4} - \\frac{1}{16}\\right)}\\\\\\frac{1}{4 \\left(- z^{2} - z - \\frac{1}{4}\\right)} & \\frac{4 z}{4 z^{2} + 4 z + 1}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[(4*z + 4)/(4*z**2 + 4*z + 1), -1/(4*(-z**2/4 - z/4 - 1/16))],\n",
       "[     1/(4*(-z**2 - z - 1/4)),        4*z/(4*z**2 + 4*z + 1)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "I want to factor this\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle z^{2} + z + \\frac{1}{4}$"
      ],
      "text/plain": [
       "z**2 + z + 1/4"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "After factoring\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\left(2 z + 1\\right)^{2}}{4}$"
      ],
      "text/plain": [
       "(2*z + 1)**2/4"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      " \n",
      " The right answer I was expecting\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(z + \\frac{1}{2}\\right)^{2}$"
      ],
      "text/plain": [
       "(z + 1/2)**2"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# v = (z**2+z)/((z-1)**2)\n",
    "# v = z/((z+1)*(z+2))\n",
    "# v = (z**2+2*z)/(z**2-3*z+2) # got the wrong answer for this one\n",
    "# # the one from the dgd\n",
    "# # v = (z*(z+1))/((z+sy.Rational(1,2))*(z+sy.Rational(1,3))**2)\n",
    "# display(v)\n",
    "# # inverse_Z_transform_partial_fraction_method(v)\n",
    "# inverse_Z_transform(v)\n",
    "\n",
    "\n",
    "## sympy factor function gives wrong answer for problems like this\n",
    "t = sy.Matrix([[z, -1], [sy.Rational(1,4), z + 1]])\n",
    "# display(t/d)\n",
    "\n",
    "display(t.inv())\n",
    "d = t.det()\n",
    "print(\"I want to factor this\")\n",
    "display(d)\n",
    "\n",
    "print(\"After factoring\")\n",
    "display(sy.factor(d))\n",
    "\n",
    "print(\" \\n The right answer I was expecting\")\n",
    "display((z+sy.Rational(1,2))**2)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#\n",
    "# State space matrix discrete time functions\n",
    "#\n",
    "def state_trans_cayley_hamilton_discrete(a):\n",
    "    '''Only implemented for a 2 by 2'''\n",
    "    a = sy.Matrix(a)\n",
    "    assert a.shape[0]==2 , \"implemented for 2by2 only\"\n",
    "    eignvals = list(a.eigenvals().keys())\n",
    "    k, a0, a1 = sy.symbols('k a0 a1')\n",
    "    eq1 = sy.Eq(eignvals[0]**k, a0+a1*eignvals[0])\n",
    "    eq2 = sy.Eq(eignvals[1]**k, a0+a1*eignvals[1])\n",
    "    display(eq1,eq2)\n",
    "    \n",
    "    e1 = eignvals[0]**k -1*(a0+a1*eignvals[0])\n",
    "    e2 = eignvals[1]**k -1*(a0+a1*eignvals[1])\n",
    "    sol = sy.solve([e1,e2], (a0,a1))\n",
    "    display(\"a0 and a1 aka alpha 0 and alpha 1\")\n",
    "    display(sy.simplify(sol[a0]),sy.simplify(sol[a1]))\n",
    "    ans = sy.simplify(sol[a0]) * sy.eye(2) + sy.simplify(sol[a1])*a\n",
    "    sy.simplify(ans)\n",
    "    display(ans)\n",
    "    \n",
    "def state_trans_matrix_discrete(a, symbol = sy.symbols('k'), show_steps = True):\n",
    "    a = sy.Matrix(a)\n",
    "    state_trans =a ** symbol\n",
    "    if show_steps:\n",
    "        display(\"the state trans matrix is\", state_trans)\n",
    "    return state_trans\n",
    "\n",
    "def state_var_forced_response_discrete(a, b,u, show_steps = True):\n",
    "    '''Takes the A matrix and the B inputs column and returns the forced \n",
    "    response of the state space.\n",
    "    '''\n",
    "    a = sy.Matrix(a)\n",
    "    b = sy.Matrix(b)\n",
    "    # T is for tau\n",
    "    k, i = sy.symbols('k i')\n",
    "    state_trans_neg = state_trans_matrix_discrete(a, k-1-i, show_steps = False)\n",
    "    state_trans = state_trans_matrix_discrete(a, k, show_steps = False)\n",
    "    # This is in case the u param is dependant on k, if it is we need to make it\n",
    "    # into an i for the summation below\n",
    "    try:\n",
    "        u = u.replace(k,i)\n",
    "    except:\n",
    "        pass\n",
    "    unevaluated_conv_sum = sy.Sum(state_trans_neg * b*u ,(i,0,k-1))\n",
    "    forced_resp = unevaluated_conv_sum.doit()\n",
    "    if show_steps:\n",
    "        display(\"Finding the forced response\")\n",
    "        display(\"A matrix\",a, \"B matrix\",b)\n",
    "        display(\"State transition matrix: \",state_trans)\n",
    "        display(\"Negative state trans matrix\", state_trans_neg)\n",
    "        display(\"The convolution summation\", unevaluated_conv_sum)\n",
    "#         display(\"Evaluated convolution summation\", conv_sum)\n",
    "        display(\"Foreced response is \", forced_resp)\n",
    "        \n",
    "    return forced_resp\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'the state trans matrix is'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}e^{a t} & 0\\\\0 & e^{b t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[exp(a*t),        0],\n",
       "[       0, exp(b*t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'the state trans matrix is'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{e^{i t w}}{2} + \\frac{e^{- i t w}}{2} & - \\frac{i e^{i t w}}{2} + \\frac{i e^{- i t w}}{2}\\\\\\frac{i e^{i t w}}{2} - \\frac{i e^{- i t w}}{2} & \\frac{e^{i t w}}{2} + \\frac{e^{- i t w}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[    exp(I*t*w)/2 + exp(-I*t*w)/2, -I*exp(I*t*w)/2 + I*exp(-I*t*w)/2],\n",
       "[I*exp(I*t*w)/2 - I*exp(-I*t*w)/2,      exp(I*t*w)/2 + exp(-I*t*w)/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note that the above is the exponent form of trig functions sine and cosine\n"
     ]
    }
   ],
   "source": [
    "a,b,w = sy.symbols('a b w')\n",
    "m = [\n",
    "     [a,0],\n",
    "    [0,b]\n",
    " ]\n",
    "state_trans_matrix(m)\n",
    "m2 = [[0,w], [-w,0]]\n",
    "state_trans_matrix(m2)\n",
    "print(\"Note that the above is the exponent form of trig functions sine and cosine\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Example of equvilant state trans matrices using different methods"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "with cayley hamilton\n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle e^{- t \\sqrt{- \\frac{g}{l}}} = a_{0} - a_{1} \\sqrt{- \\frac{g}{l}}$"
      ],
      "text/plain": [
       "Eq(exp(-t*sqrt(-g/l)), a0 - a1*sqrt(-g/l))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle e^{t \\sqrt{- \\frac{g}{l}}} = a_{0} + a_{1} \\sqrt{- \\frac{g}{l}}$"
      ],
      "text/plain": [
       "Eq(exp(t*sqrt(-g/l)), a0 + a1*sqrt(-g/l))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'a0 and a1 aka alpha 0 and alpha 1'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\cosh{\\left(t \\sqrt{- \\frac{g}{l}} \\right)}$"
      ],
      "text/plain": [
       "cosh(t*sqrt(-g/l))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{\\sinh{\\left(t \\sqrt{- \\frac{g}{l}} \\right)}}{\\sqrt{- \\frac{g}{l}}}$"
      ],
      "text/plain": [
       "sinh(t*sqrt(-g/l))/sqrt(-g/l)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\cosh{\\left(t \\sqrt{- \\frac{g}{l}} \\right)} & - \\frac{g \\sinh{\\left(t \\sqrt{- \\frac{g}{l}} \\right)}}{l \\sqrt{- \\frac{g}{l}}}\\\\\\frac{\\sinh{\\left(t \\sqrt{- \\frac{g}{l}} \\right)}}{\\sqrt{- \\frac{g}{l}}} & \\cosh{\\left(t \\sqrt{- \\frac{g}{l}} \\right)}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[           cosh(t*sqrt(-g/l)), -g*sinh(t*sqrt(-g/l))/(l*sqrt(-g/l))],\n",
       "[sinh(t*sqrt(-g/l))/sqrt(-g/l),                   cosh(t*sqrt(-g/l))]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'We can also do this for an equivalent answer'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{e^{t \\sqrt{- \\frac{g}{l}}}}{2} + \\frac{e^{- t \\sqrt{- \\frac{g}{l}}}}{2} & - \\frac{g e^{t \\sqrt{- \\frac{g}{l}}}}{2 l \\sqrt{- \\frac{g}{l}}} + \\frac{g e^{- t \\sqrt{- \\frac{g}{l}}}}{2 l \\sqrt{- \\frac{g}{l}}}\\\\- \\frac{l \\sqrt{- \\frac{g}{l}} e^{t \\sqrt{- \\frac{g}{l}}}}{2 g} + \\frac{l \\sqrt{- \\frac{g}{l}} e^{- t \\sqrt{- \\frac{g}{l}}}}{2 g} & \\frac{e^{t \\sqrt{- \\frac{g}{l}}}}{2} + \\frac{e^{- t \\sqrt{- \\frac{g}{l}}}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                   exp(t*sqrt(-g/l))/2 + exp(-t*sqrt(-g/l))/2, -g*exp(t*sqrt(-g/l))/(2*l*sqrt(-g/l)) + g*exp(-t*sqrt(-g/l))/(2*l*sqrt(-g/l))],\n",
       "[-l*sqrt(-g/l)*exp(t*sqrt(-g/l))/(2*g) + l*sqrt(-g/l)*exp(-t*sqrt(-g/l))/(2*g),                                    exp(t*sqrt(-g/l))/2 + exp(-t*sqrt(-g/l))/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'the state trans matrix is'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{e^{t \\sqrt{- \\frac{g}{l}}}}{2} + \\frac{e^{- t \\sqrt{- \\frac{g}{l}}}}{2} & - \\frac{g e^{t \\sqrt{- \\frac{g}{l}}}}{2 l \\sqrt{- \\frac{g}{l}}} + \\frac{g e^{- t \\sqrt{- \\frac{g}{l}}}}{2 l \\sqrt{- \\frac{g}{l}}}\\\\- \\frac{l \\sqrt{- \\frac{g}{l}} e^{t \\sqrt{- \\frac{g}{l}}}}{2 g} + \\frac{l \\sqrt{- \\frac{g}{l}} e^{- t \\sqrt{- \\frac{g}{l}}}}{2 g} & \\frac{e^{t \\sqrt{- \\frac{g}{l}}}}{2} + \\frac{e^{- t \\sqrt{- \\frac{g}{l}}}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                   exp(t*sqrt(-g/l))/2 + exp(-t*sqrt(-g/l))/2, -g*exp(t*sqrt(-g/l))/(2*l*sqrt(-g/l)) + g*exp(-t*sqrt(-g/l))/(2*l*sqrt(-g/l))],\n",
       "[-l*sqrt(-g/l)*exp(t*sqrt(-g/l))/(2*g) + l*sqrt(-g/l)*exp(-t*sqrt(-g/l))/(2*g),                                    exp(t*sqrt(-g/l))/2 + exp(-t*sqrt(-g/l))/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Finding the homogenous response'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'A matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & - \\frac{g}{l}\\\\1 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0, -g/l],\n",
       "[1,    0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Initial condition'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'x homogenous response at time t is '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{g e^{t \\sqrt{- \\frac{g}{l}}}}{2 l \\sqrt{- \\frac{g}{l}}} + \\frac{g e^{- t \\sqrt{- \\frac{g}{l}}}}{2 l \\sqrt{- \\frac{g}{l}}}\\\\\\frac{e^{t \\sqrt{- \\frac{g}{l}}}}{2} + \\frac{e^{- t \\sqrt{- \\frac{g}{l}}}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-g*exp(t*sqrt(-g/l))/(2*l*sqrt(-g/l)) + g*exp(-t*sqrt(-g/l))/(2*l*sqrt(-g/l))],\n",
       "[                                   exp(t*sqrt(-g/l))/2 + exp(-t*sqrt(-g/l))/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Note that cosh(imaginary * x) =  cos (x), since the above has sqrt(-1) we can also write them in terms of sines and cosines and remove the imaginary term. We can see from below that their series equivlant is the same \n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "'These two are equal, We can also see that their series below are equal'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\cos{\\left(t \\sqrt{\\frac{g}{l}} \\right)}$"
      ],
      "text/plain": [
       "cos(t*sqrt(g/l))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\cosh{\\left(t \\sqrt{- \\frac{g}{l}} \\right)}$"
      ],
      "text/plain": [
       "cosh(t*sqrt(-g/l))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1 - \\frac{g t^{2}}{2 l} + O\\left(t^{3}\\right)$"
      ],
      "text/plain": [
       "1 - g*t**2/(2*l) + O(t**3)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle 1 - \\frac{g t^{2}}{2 l} + O\\left(t^{3}\\right)$"
      ],
      "text/plain": [
       "1 - g*t**2/(2*l) + O(t**3)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "g,l,t = sy.symbols('g l t')\n",
    "a = sy.Matrix([\n",
    "    [0,-g/l],\n",
    "    [1,0]\n",
    "])\n",
    "# a = sy.Matrix([\n",
    "#     [0,1],[-2,-3]\n",
    "# ])\n",
    "# state_trans_matrix_inverse_laplace_method(a)\n",
    "print(\"with cayley hamilton\")\n",
    "state_trans_cayley_hamilton_method(a)\n",
    "display(\"We can also do this for an equivalent answer\", sy.exp(a*t))\n",
    "state_var_homo_response(a,[0,1])\n",
    "print(\"Note that cosh(imaginary * x) =  cos (x), since the above has sqrt(-1) we can also write them in terms of sines and cosines and remove the imaginary term. We can see from below that their series equivlant is the same \")\n",
    "cos_eq = sy.cos(sy.sqrt(g/l)*t)\n",
    "cosh_eq = sy.cosh(sy.sqrt(-g/l)*t)\n",
    "display(\"These two are equal, We can also see that their series below are equal\", cos_eq, cosh_eq)\n",
    "display(sy.series(cos_eq,t,n=3))\n",
    "display(sy.series(cosh_eq,t,n=3))\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'finding the output response'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'A matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-2 & 0\\\\1 & -1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-2,  0],\n",
       "[ 1, -1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'B matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'C matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}2 & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[2, 1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'D matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'initial condition'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'homogenous resp'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'forced resp'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) e^{- 2 t}\\\\\\left(e^{- t} - e^{- 2 t}\\right) \\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) + \\left(- \\frac{5 e^{2 t}}{2} + 5 e^{t} - \\frac{5}{2}\\right) e^{- t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                                       (5*exp(2*t)/2 - 5/2)*exp(-2*t)],\n",
       "[(exp(-t) - exp(-2*t))*(5*exp(2*t)/2 - 5/2) + (-5*exp(2*t)/2 + 5*exp(t) - 5/2)*exp(-t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Full x response ie x homogenous + x forced'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) e^{- 2 t}\\\\\\left(e^{- t} - e^{- 2 t}\\right) \\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) + \\left(- \\frac{5 e^{2 t}}{2} + 5 e^{t} - \\frac{5}{2}\\right) e^{- t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                                                       (5*exp(2*t)/2 - 5/2)*exp(-2*t)],\n",
       "[(exp(-t) - exp(-2*t))*(5*exp(2*t)/2 - 5/2) + (-5*exp(2*t)/2 + 5*exp(t) - 5/2)*exp(-t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Full response times c col at time t (or k if discrete)'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\left(e^{- t} - e^{- 2 t}\\right) \\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) + 2 \\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) e^{- 2 t} + \\left(- \\frac{5 e^{2 t}}{2} + 5 e^{t} - \\frac{5}{2}\\right) e^{- t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[(exp(-t) - exp(-2*t))*(5*exp(2*t)/2 - 5/2) + 2*(5*exp(2*t)/2 - 5/2)*exp(-2*t) + (-5*exp(2*t)/2 + 5*exp(t) - 5/2)*exp(-t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'the output y is '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\left(e^{- t} - e^{- 2 t}\\right) \\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) + 2 \\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) e^{- 2 t} + \\left(- \\frac{5 e^{2 t}}{2} + 5 e^{t} - \\frac{5}{2}\\right) e^{- t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[(exp(-t) - exp(-2*t))*(5*exp(2*t)/2 - 5/2) + 2*(5*exp(2*t)/2 - 5/2)*exp(-2*t) + (-5*exp(2*t)/2 + 5*exp(t) - 5/2)*exp(-t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\left(e^{- t} - e^{- 2 t}\\right) \\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) + 2 \\left(\\frac{5 e^{2 t}}{2} - \\frac{5}{2}\\right) e^{- 2 t} + \\left(- \\frac{5 e^{2 t}}{2} + 5 e^{t} - \\frac{5}{2}\\right) e^{- t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[(exp(-t) - exp(-2*t))*(5*exp(2*t)/2 - 5/2) + 2*(5*exp(2*t)/2 - 5/2)*exp(-2*t) + (-5*exp(2*t)/2 + 5*exp(t) - 5/2)*exp(-t)]])"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# example of a state space model solution\n",
    "# x' = Ax + Bu\n",
    "# y  = Cx + Du\n",
    "\n",
    "a= sy.Matrix([\n",
    "    [-2,0],\n",
    "    [1,-1]\n",
    "])\n",
    "# ic = sy.Matrix([[2],[3]])\n",
    "ic = [0,0]\n",
    "b = sy.Matrix([\n",
    "    [1],\n",
    "    [0]\n",
    "])\n",
    "u = 5\n",
    "\n",
    "c = [\n",
    "    [2,1]\n",
    "]\n",
    "\n",
    "# state_var_homo_response(a, ic)\n",
    "\n",
    "output_response(a,b,c,[0],u, ic)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Finding the forced response'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'A matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{3}{4} & - \\frac{1}{4}\\\\- \\frac{1}{4} & - \\frac{3}{4}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-3/4, -1/4],\n",
       "[-1/4, -3/4]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'B matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'State transition matrix: '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\left(-1\\right)^{k}}{2} + \\frac{\\left(- \\frac{1}{2}\\right)^{k}}{2} & \\frac{\\left(-1\\right)^{k}}{2} - \\frac{\\left(- \\frac{1}{2}\\right)^{k}}{2}\\\\\\frac{\\left(-1\\right)^{k}}{2} - \\frac{\\left(- \\frac{1}{2}\\right)^{k}}{2} & \\frac{\\left(-1\\right)^{k}}{2} + \\frac{\\left(- \\frac{1}{2}\\right)^{k}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[(-1)**k/2 + (-1/2)**k/2, (-1)**k/2 - (-1/2)**k/2],\n",
       "[(-1)**k/2 - (-1/2)**k/2, (-1)**k/2 + (-1/2)**k/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Negative state trans matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{\\left(-1\\right)^{- i + k - 1}}{2} + \\frac{\\left(- \\frac{1}{2}\\right)^{- i + k - 1}}{2} & \\frac{\\left(-1\\right)^{- i + k - 1}}{2} - \\frac{\\left(- \\frac{1}{2}\\right)^{- i + k - 1}}{2}\\\\\\frac{\\left(-1\\right)^{- i + k - 1}}{2} - \\frac{\\left(- \\frac{1}{2}\\right)^{- i + k - 1}}{2} & \\frac{\\left(-1\\right)^{- i + k - 1}}{2} + \\frac{\\left(- \\frac{1}{2}\\right)^{- i + k - 1}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[(-1)**(-i + k - 1)/2 + (-1/2)**(-i + k - 1)/2, (-1)**(-i + k - 1)/2 - (-1/2)**(-i + k - 1)/2],\n",
       "[(-1)**(-i + k - 1)/2 - (-1/2)**(-i + k - 1)/2, (-1)**(-i + k - 1)/2 + (-1/2)**(-i + k - 1)/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'The convolution summation'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sum_{i=0}^{k - 1} \\left[\\begin{matrix}\\frac{\\left(-1\\right)^{- i + k - 1}}{2} + \\frac{\\left(- \\frac{1}{2}\\right)^{- i + k - 1}}{2}\\\\\\frac{\\left(-1\\right)^{- i + k - 1}}{2} - \\frac{\\left(- \\frac{1}{2}\\right)^{- i + k - 1}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Sum(Matrix([\n",
       "[(-1)**(-i + k - 1)/2 + (-1/2)**(-i + k - 1)/2],\n",
       "[(-1)**(-i + k - 1)/2 - (-1/2)**(-i + k - 1)/2]]), (i, 0, k - 1))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Foreced response is '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(-1\\right)^{k} \\left(\\frac{1}{2} - \\frac{\\left(-1\\right)^{k}}{2}\\right)}{2} - \\left(- \\frac{1}{2}\\right)^{k} \\left(\\frac{1}{3} - \\frac{\\left(-2\\right)^{k}}{3}\\right)\\\\- \\frac{\\left(-1\\right)^{k} \\left(\\frac{1}{2} - \\frac{\\left(-1\\right)^{k}}{2}\\right)}{2} + \\left(- \\frac{1}{2}\\right)^{k} \\left(\\frac{1}{3} - \\frac{\\left(-2\\right)^{k}}{3}\\right)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-(-1)**k*(1/2 - (-1)**k/2)/2 - (-1/2)**k*(1/3 - (-2)**k/3)],\n",
       "[-(-1)**k*(1/2 - (-1)**k/2)/2 + (-1/2)**k*(1/3 - (-2)**k/3)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'finding the output response'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'A matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{3}{4} & - \\frac{1}{4}\\\\- \\frac{1}{4} & - \\frac{3}{4}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-3/4, -1/4],\n",
       "[-1/4, -3/4]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'B matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'C matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[1, 0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'D matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}4\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[4]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'initial condition'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'homogenous resp'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'forced resp'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(-1\\right)^{k} \\left(\\frac{1}{2} - \\frac{\\left(-1\\right)^{k}}{2}\\right)}{2} - \\left(- \\frac{1}{2}\\right)^{k} \\left(\\frac{1}{3} - \\frac{\\left(-2\\right)^{k}}{3}\\right)\\\\- \\frac{\\left(-1\\right)^{k} \\left(\\frac{1}{2} - \\frac{\\left(-1\\right)^{k}}{2}\\right)}{2} + \\left(- \\frac{1}{2}\\right)^{k} \\left(\\frac{1}{3} - \\frac{\\left(-2\\right)^{k}}{3}\\right)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-(-1)**k*(1/2 - (-1)**k/2)/2 - (-1/2)**k*(1/3 - (-2)**k/3)],\n",
       "[-(-1)**k*(1/2 - (-1)**k/2)/2 + (-1/2)**k*(1/3 - (-2)**k/3)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Full x response ie x homogenous + x forced'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(-1\\right)^{k} \\left(\\frac{1}{2} - \\frac{\\left(-1\\right)^{k}}{2}\\right)}{2} - \\left(- \\frac{1}{2}\\right)^{k} \\left(\\frac{1}{3} - \\frac{\\left(-2\\right)^{k}}{3}\\right)\\\\- \\frac{\\left(-1\\right)^{k} \\left(\\frac{1}{2} - \\frac{\\left(-1\\right)^{k}}{2}\\right)}{2} + \\left(- \\frac{1}{2}\\right)^{k} \\left(\\frac{1}{3} - \\frac{\\left(-2\\right)^{k}}{3}\\right)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-(-1)**k*(1/2 - (-1)**k/2)/2 - (-1/2)**k*(1/3 - (-2)**k/3)],\n",
       "[-(-1)**k*(1/2 - (-1)**k/2)/2 + (-1/2)**k*(1/3 - (-2)**k/3)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Full response times c col at time t (or k if discrete)'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(-1\\right)^{k} \\left(\\frac{1}{2} - \\frac{\\left(-1\\right)^{k}}{2}\\right)}{2} - \\left(- \\frac{1}{2}\\right)^{k} \\left(\\frac{1}{3} - \\frac{\\left(-2\\right)^{k}}{3}\\right)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-(-1)**k*(1/2 - (-1)**k/2)/2 - (-1/2)**k*(1/3 - (-2)**k/3)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'the output y is '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{\\left(-1\\right)^{k} \\left(\\frac{1}{2} - \\frac{\\left(-1\\right)^{k}}{2}\\right)}{2} - \\left(- \\frac{1}{2}\\right)^{k} \\left(\\frac{1}{3} - \\frac{\\left(-2\\right)^{k}}{3}\\right) + 4\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-(-1)**k*(1/2 - (-1)**k/2)/2 - (-1/2)**k*(1/3 - (-2)**k/3) + 4]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "# Obtaining the state transition matrix using different methods yields\n",
    "# equivelant results\n",
    "a = sy.Matrix([\n",
    "[sy.Rational(-3,4) , sy.Rational(-1,4)],\n",
    "    [sy.Rational(-1,4) , sy.Rational(-3,4)]\n",
    "])\n",
    "b = [\n",
    "    [1],\n",
    "    [0]\n",
    "]\n",
    "u=1\n",
    "c = [[1,0]]\n",
    "d = [4]\n",
    "ic = [[0],\n",
    "     [0]]\n",
    "# matrix_info(a)\n",
    "# state_trans = state_trans_matrix(a)\n",
    "# state_trans_cayley_hamilton_method(a)\n",
    "state_var_forced_response_discrete(a, b, u)\n",
    "# state_var_homo_response(a, ic)\n",
    "my_ans = sy.simplify(output_response(a,b,c,d,u,ic, discrete = True))\n",
    "# display(state_trans_discrete_cayley_hamilton(a))\n",
    "# state_trans_matrix_discrete(a)\n",
    "\n",
    "######\n",
    "# # comparing my solution against the solution manual, and mine is right!\n",
    "# k = sy.symbols('k')\n",
    "# sy.collect(sy.expand(sy.simplify(a**k)[0]),k).equals((-1)**k * (2**k +1)/2**(k+1))\n",
    "# display(sy.expand(sy.simplify(a**k)[0]).equals((-1)**k * (2**k +1)/2**(k+1)))\n",
    "\n",
    "# i = sy.symbols('i')\n",
    "# book_ans = 4 + sy.Sum((-1)**(k-1-i)*(2**(k-1-i)+1)/(2**(k-i)), (i,0,k-1))\n",
    "# display(\"answer in the book\",book_ans)\n",
    "# display(\"answer in the book\",book_ans)\n",
    "# book_ans = book_ans.doit()\n",
    "# display(\"answer in the book\",book_ans)\n",
    "\n",
    "# display(\"my answer\",my_ans[0])\n",
    "# # display(book_ans.equals(my_ans[0]))\n",
    "# print(my_ans[0].subs(k,1).equals(book_ans.subs(k,1)))\n",
    "######"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "A-Lamda x I: \n"
     ]
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- L & 0\\\\t & - L\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-L,  0],\n",
       "[ t, -L]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "{'det': 0, 'inverse': 'Not invertable, det=0', 'charFunction': 'L^{2}'}"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Sympy eigenvalues: '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[0]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Sympy eigenvectors: '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "[[0, 1]]"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'the state trans matrix is'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0\\\\t^{2} & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[   1, 0],\n",
       "[t**2, 1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Finding the forced response'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'A matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0\\\\t & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0, 0],\n",
       "[t, 0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'B matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1\\\\t\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1],\n",
       "[t]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'State transition matrix: '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0\\\\t^{2} & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[   1, 0],\n",
       "[t**2, 1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Negative state trans matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0\\\\- T t & 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[   1, 0],\n",
       "[-T*t, 1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'The convolution integral'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\int\\limits_{0}^{t} \\left[\\begin{matrix}T\\\\T \\left(- T t + t\\right)\\end{matrix}\\right]\\, dT$"
      ],
      "text/plain": [
       "Integral(Matrix([\n",
       "[           T],\n",
       "[T*(-T*t + t)]]), (T, 0, t))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Evaluated convolution integral'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{t^{2}}{2}\\\\- \\frac{t^{4}}{3} + \\frac{t^{3}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[          t**2/2],\n",
       "[-t**4/3 + t**3/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Foreced response is '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{t^{2}}{2}\\\\\\frac{t^{4}}{6} + \\frac{t^{3}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[         t**2/2],\n",
       "[t**4/6 + t**3/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'finding the output response'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'A matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 0\\\\t & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0, 0],\n",
       "[t, 0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'B matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1\\\\t\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1],\n",
       "[t]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'C matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & t\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[1, t]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'D matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}t\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[t]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'initial condition'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1\\\\1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1],\n",
       "[1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'homogenous resp'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1\\\\t^{2} + 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[       1],\n",
       "[t**2 + 1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'forced resp'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{t^{2}}{2}\\\\\\frac{t^{4}}{6} + \\frac{t^{3}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[         t**2/2],\n",
       "[t**4/6 + t**3/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Full x response ie x homogenous + x forced'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{t^{2}}{2} + 1\\\\\\frac{t^{4}}{6} + \\frac{t^{3}}{2} + t^{2} + 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[                t**2/2 + 1],\n",
       "[t**4/6 + t**3/2 + t**2 + 1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Full response times c col at time t (or k if discrete)'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{t^{2}}{2} + t \\left(t^{2} + 1\\right) + t \\left(\\frac{t^{4}}{6} + \\frac{t^{3}}{2}\\right) + 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[t**2/2 + t*(t**2 + 1) + t*(t**4/6 + t**3/2) + 1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'the output y is '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{3 t^{2}}{2} + t \\left(t^{2} + 1\\right) + t \\left(\\frac{t^{4}}{6} + \\frac{t^{3}}{2}\\right) + 1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[3*t**2/2 + t*(t**2 + 1) + t*(t**4/6 + t**3/2) + 1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "### Q5\n",
    "t, T = sy.symbols('t T')\n",
    "a = [\n",
    "    [0,0],\n",
    "    [t,0]\n",
    "]\n",
    "b =  sy.Matrix([\n",
    "    [1],\n",
    "    [t]\n",
    "])\n",
    "u = t\n",
    "c = [[1,t]]\n",
    "d=[t]\n",
    "ic = [[1],[1]]\n",
    "\n",
    "matrix_info(a)\n",
    "get_eign(a)\n",
    "state_trans = state_trans_matrix(a)\n",
    "state_var_forced_response(a, b,u)\n",
    "\n",
    "# state_var_homo_response(a, ic)\n",
    "my_ans = sy.simplify(output_response(a,b,c,d,u,ic))\n",
    "\n",
    "### Checking my answer against the solution manual\n",
    "# wtf the sol manual gave us two solutions to the \n",
    "# problem but they're both WRONG and not equal to each other\n",
    "\n",
    "# book_sol =  1+t+3*t**2/2 + t**3/2 + t**4/3 + t**5/8\n",
    "# book_sol2 = 3*t**2/2 + t**5/8 + t**4/3\n",
    "# print(my_ans[0].equals(book_sol))\n",
    "# print(my_ans[0].equals(book_sol2))\n",
    "# book_sol.subs(t,1).equals(book_sol2.subs(t,1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Example usage"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'For '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1 & 0 & 0\\\\0 & -4 & 4\\\\0 & -1 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-1,  0, 0],\n",
       "[ 0, -4, 4],\n",
       "[ 0, -1, 0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'the state trans matrix is'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}e^{- t} & 0 & 0\\\\0 & - 2 t e^{- 2 t} + e^{- 2 t} & 4 t e^{- 2 t}\\\\0 & - t e^{- 2 t} & 2 t e^{- 2 t} + e^{- 2 t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[exp(-t),                          0,                         0],\n",
       "[      0, -2*t*exp(-2*t) + exp(-2*t),             4*t*exp(-2*t)],\n",
       "[      0,               -t*exp(-2*t), 2*t*exp(-2*t) + exp(-2*t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'For '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}x & y\\\\- y & x\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[ x, y],\n",
       "[-y, x]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'the state trans matrix is'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\frac{e^{t x - i t y}}{2} + \\frac{e^{t x + i t y}}{2} & \\frac{i e^{t x - i t y}}{2} - \\frac{i e^{t x + i t y}}{2}\\\\- \\frac{i e^{t x - i t y}}{2} + \\frac{i e^{t x + i t y}}{2} & \\frac{e^{t x - i t y}}{2} + \\frac{e^{t x + i t y}}{2}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[     exp(t*x - I*t*y)/2 + exp(t*x + I*t*y)/2, I*exp(t*x - I*t*y)/2 - I*exp(t*x + I*t*y)/2],\n",
       "[-I*exp(t*x - I*t*y)/2 + I*exp(t*x + I*t*y)/2,     exp(t*x - I*t*y)/2 + exp(t*x + I*t*y)/2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "# finding state transition matrix\n",
    "a = sy.Matrix([\n",
    "    [-1,0,0],\n",
    "    [0,-4,4],\n",
    "    [0,-1,0]\n",
    "])\n",
    "display(\"For \",a)\n",
    "ans = state_trans_matrix(a)\n",
    "x,y = sy.symbols('x y')\n",
    "a = sy.Matrix([\n",
    "    [x,y],\n",
    "    [-y,x]\n",
    "])\n",
    "display(\"For \",a)\n",
    "\n",
    "ans = state_trans_matrix(a)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finding the output for the following:\n",
    "\n",
    "\n",
    "$$\\dot{x} = \\begin{bmatrix} -2 & 0 \\\\ 0 & -3 \\end{bmatrix}\n",
    "\\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix}\n",
    "+\n",
    "\\begin{bmatrix} 1 \\\\ 1 \\end{bmatrix}\n",
    "u\n",
    "\\\\\n",
    "y = [-1,2]\\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix}\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'finding the output response'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'A matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-2 & 0\\\\0 & -3\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-2,  0],\n",
       "[ 0, -3]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'B matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1\\\\1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[1],\n",
       "[1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'C matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}-1 & 2\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-1, 2]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'D matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'initial condition'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'homogenous resp'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'forced resp'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\left(\\frac{e^{2 t}}{2} - \\frac{1}{2}\\right) e^{- 2 t}\\\\\\left(\\frac{e^{3 t}}{3} - \\frac{1}{3}\\right) e^{- 3 t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[(exp(2*t)/2 - 1/2)*exp(-2*t)],\n",
       "[(exp(3*t)/3 - 1/3)*exp(-3*t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Full x response ie x homogenous + x forced'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}\\left(\\frac{e^{2 t}}{2} - \\frac{1}{2}\\right) e^{- 2 t}\\\\\\left(\\frac{e^{3 t}}{3} - \\frac{1}{3}\\right) e^{- 3 t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[(exp(2*t)/2 - 1/2)*exp(-2*t)],\n",
       "[(exp(3*t)/3 - 1/3)*exp(-3*t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Full response times c col at time t (or k if discrete)'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\left(\\frac{e^{2 t}}{2} - \\frac{1}{2}\\right) e^{- 2 t} + 2 \\left(\\frac{e^{3 t}}{3} - \\frac{1}{3}\\right) e^{- 3 t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-(exp(2*t)/2 - 1/2)*exp(-2*t) + 2*(exp(3*t)/3 - 1/3)*exp(-3*t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'the output y is '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\left(\\frac{e^{2 t}}{2} - \\frac{1}{2}\\right) e^{- 2 t} + 2 \\left(\\frac{e^{3 t}}{3} - \\frac{1}{3}\\right) e^{- 3 t}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([[-(exp(2*t)/2 - 1/2)*exp(-2*t) + 2*(exp(3*t)/3 - 1/3)*exp(-3*t)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\frac{1}{6} + \\frac{e^{- 2 t}}{2} - \\frac{2 e^{- 3 t}}{3}$"
      ],
      "text/plain": [
       "1/6 + exp(-2*t)/2 - 2*exp(-3*t)/3"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "a = [\n",
    "    [-2,0],\n",
    "    [0,-3]\n",
    "]\n",
    "b = [[1],[1]]\n",
    "c = [[-1,2]]\n",
    "d = [0]\n",
    "u=1\n",
    "# initial condition\n",
    "ic = [[0],[0]]\n",
    "sy.expand(output_response(a,b,c,d,u,ic)[0])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Discrete time example\n",
    "\n",
    "\n",
    "$$ X(k+1)  =  A(k)X(k) + B(k)U(k)\\\\\n",
    "Y(k) = C(k)X(k) + D(k)U(k)\n",
    "$$\n",
    "\n",
    "<br>\n",
    "\n",
    "---\n",
    "\n",
    "find the output for the following when k>=0 and u(k) = 1\n",
    "\n",
    "$$ X(k) = \\begin{bmatrix} 1 & 0 \\\\ k+1 & 0 \\end{bmatrix}\n",
    "\\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix}\n",
    "+\n",
    "\\begin{bmatrix} 0 \\\\ 0 \\end{bmatrix}\n",
    "[u(k)]\n",
    "\\\\\n",
    "y = [1,\\frac{1}{k+1}]\\begin{bmatrix} x_1 \\\\ x_2 \\end{bmatrix}\n",
    "+[1][u(k)]\n",
    "\\\\\n",
    "X(0) = \\begin{bmatrix} 1 \\\\ 1 \\end{bmatrix}\n",
    "$$\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Finding the forced response'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'A matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1 & 0\\\\k + 1 & 0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[    1, 0],\n",
       "[k + 1, 0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'B matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[0]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'State transition matrix: '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(\\left[\\begin{matrix}1 & 0\\\\k + 1 & 0\\end{matrix}\\right]\\right)^{k}$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[    1, 0],\n",
       "[k + 1, 0]])**k"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Negative state trans matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left(\\left[\\begin{matrix}1 & 0\\\\k + 1 & 0\\end{matrix}\\right]\\right)^{- i + k - 1}$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[    1, 0],\n",
       "[k + 1, 0]])**(-i + k - 1)"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'The convolution summation'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sum_{i=0}^{k - 1} \\left(\\left[\\begin{matrix}1 & 0\\\\k + 1 & 0\\end{matrix}\\right]\\right)^{- i + k - 1} \\left[\\begin{matrix}0\\\\0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Sum(Matrix([\n",
       "[    1, 0],\n",
       "[k + 1, 0]])**(-i + k - 1)*Matrix([\n",
       "[0],\n",
       "[0]]), (i, 0, k - 1))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Foreced response is '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "None"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "\n",
    "k = sy.symbols('k')\n",
    "a = [\n",
    "    [1 , 0 ],\n",
    "    [k+1, 0]\n",
    "]\n",
    "b = [\n",
    "    [0],\n",
    "    [0]\n",
    "    ]\n",
    "c = [\n",
    "    [1, 1/(k+1)]\n",
    "]\n",
    "d = [1]\n",
    "u = 1\n",
    "ic = [\n",
    "    [1],\n",
    "    [1]\n",
    "]\n",
    "# output_response(a,b,c,d,u,ic, discrete=True)\n",
    "state_var_forced_response_discrete(a,b,u)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Finding the forced response'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'A matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0 & 1\\\\-0.21 & -1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[    0,  1],\n",
       "[-0.21, -1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'B matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\1\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[0],\n",
       "[1]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'State transition matrix: '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1.75 \\left(-0.3\\right)^{k} - 0.75 \\left(-0.7\\right)^{k} & 2.5 \\left(-0.3\\right)^{k} - 2.5 \\left(-0.7\\right)^{k}\\\\- 0.525 \\left(-0.3\\right)^{k} + 0.525 \\left(-0.7\\right)^{k} & - 0.75 \\left(-0.3\\right)^{k} + 1.75 \\left(-0.7\\right)^{k}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[   1.75*(-0.3)**k - 0.75*(-0.7)**k,    2.5*(-0.3)**k - 2.5*(-0.7)**k],\n",
       "[-0.525*(-0.3)**k + 0.525*(-0.7)**k, -0.75*(-0.3)**k + 1.75*(-0.7)**k]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Negative state trans matrix'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}1.75 \\left(-0.3\\right)^{- i + k - 1} - 0.75 \\left(-0.7\\right)^{- i + k - 1} & 2.5 \\left(-0.3\\right)^{- i + k - 1} - 2.5 \\left(-0.7\\right)^{- i + k - 1}\\\\- 0.525 \\left(-0.3\\right)^{- i + k - 1} + 0.525 \\left(-0.7\\right)^{- i + k - 1} & - 0.75 \\left(-0.3\\right)^{- i + k - 1} + 1.75 \\left(-0.7\\right)^{- i + k - 1}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[   1.75*(-0.3)**(-i + k - 1) - 0.75*(-0.7)**(-i + k - 1),    2.5*(-0.3)**(-i + k - 1) - 2.5*(-0.7)**(-i + k - 1)],\n",
       "[-0.525*(-0.3)**(-i + k - 1) + 0.525*(-0.7)**(-i + k - 1), -0.75*(-0.3)**(-i + k - 1) + 1.75*(-0.7)**(-i + k - 1)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'The convolution summation'"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\sum_{i=0}^{k - 1} \\left[\\begin{matrix}\\left(-1\\right)^{i} \\left(2.5 \\left(-0.3\\right)^{- i + k - 1} - 2.5 \\left(-0.7\\right)^{- i + k - 1}\\right)\\\\\\left(-1\\right)^{i} \\left(- 0.75 \\left(-0.3\\right)^{- i + k - 1} + 1.75 \\left(-0.7\\right)^{- i + k - 1}\\right)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Sum(Matrix([\n",
       "[   (-1)**i*(2.5*(-0.3)**(-i + k - 1) - 2.5*(-0.7)**(-i + k - 1))],\n",
       "[(-1)**i*(-0.75*(-0.3)**(-i + k - 1) + 1.75*(-0.7)**(-i + k - 1))]]), (i, 0, k - 1))"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": [
       "'Foreced response is '"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- 3.57142857142857 \\left(-0.3\\right)^{k} \\left(-1.0 + 1.0 \\left(-0.3\\right)^{- k} \\left(-1\\right)^{k}\\right) + 8.33333333333333 \\left(-0.7\\right)^{k} \\left(-1.0 + 1.0 \\left(-0.7\\right)^{- k} \\left(-1\\right)^{k}\\right)\\\\1.07142857142857 \\left(-0.3\\right)^{k} \\left(-1.0 + 1.0 \\left(-0.3\\right)^{- k} \\left(-1\\right)^{k}\\right) - 5.83333333333333 \\left(-0.7\\right)^{k} \\left(-1.0 + 1.0 \\left(-0.7\\right)^{- k} \\left(-1\\right)^{k}\\right)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-3.57142857142857*(-0.3)**k*(-1.0 + 1.0*(-0.3)**(-k)*(-1)**k) + 8.33333333333333*(-0.7)**k*(-1.0 + 1.0*(-0.7)**(-k)*(-1)**k)],\n",
       "[ 1.07142857142857*(-0.3)**k*(-1.0 + 1.0*(-0.3)**(-k)*(-1)**k) - 5.83333333333333*(-0.7)**k*(-1.0 + 1.0*(-0.7)**(-k)*(-1)**k)]])"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- 3.57 \\left(-0.3\\right)^{k} \\left(-1.0 + 1.0 \\left(-0.3\\right)^{- k} \\left(-1\\right)^{k}\\right) + 8.33 \\left(-0.7\\right)^{k} \\left(-1.0 + 1.0 \\left(-0.7\\right)^{- k} \\left(-1\\right)^{k}\\right)\\\\1.07 \\left(-0.3\\right)^{k} \\left(-1.0 + 1.0 \\left(-0.3\\right)^{- k} \\left(-1\\right)^{k}\\right) - 5.83 \\left(-0.7\\right)^{k} \\left(-1.0 + 1.0 \\left(-0.7\\right)^{- k} \\left(-1\\right)^{k}\\right)\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-3.57*(-0.3)**k*(-1.0 + 1.0*(-0.3)**(-k)*(-1)**k) + 8.33*(-0.7)**k*(-1.0 + 1.0*(-0.7)**(-k)*(-1)**k)],\n",
       "[ 1.07*(-0.3)**k*(-1.0 + 1.0*(-0.3)**(-k)*(-1)**k) - 5.83*(-0.7)**k*(-1.0 + 1.0*(-0.7)**(-k)*(-1)**k)]])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# example for solving discrete SS\n",
    "# pdf link: https://nptel.ac.in/content/storage2/courses/108103008/PDF/module7/m7_lec4.pdf\n",
    "a = sy.Matrix([\n",
    "    [0, 1],\n",
    "    [-0.21, -1]\n",
    "])\n",
    "b = [[0],[1]]\n",
    "u = (-1)**k\n",
    "d = [0]\n",
    "c = [[0,1]]\n",
    "ic = [[1],[0]]\n",
    "# display(state_trans_discrete_cayley_hamilton(a))\n",
    "# state_trans_discrete_matrix(a)\n",
    "round_expr(state_var_forced_response_discrete(a,b,u),2)\n",
    "\n",
    "# my_y = sy.expand(round_expr(output_response(a,b,c,d,u,ic,discrete=True, show_steps=True),2)[0], k)\n",
    "# y_from_pdf = 0.475*(-3)**k - 5.3*(-0.7)**k + (-0.3)**k*(3.33)**k+5.825*(-0.7)**(k-1)*(1.43)**k\n",
    "# display(y_from_pdf)\n",
    "\n",
    "# display(\"If we compare the output from the pdf answer and my answer when we set k = 1 we get\")\n",
    "# display(\"pdf answer\",y_from_pdf.subs(k,1))\n",
    "# display(\"my answer\",my_y.subs(k,1))\n",
    "# display(\" I think that my answer is right and their answer is wrong, because my function solved other questions correctly\")"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "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.8.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}