docstr: Calculate time-domain solution for
  de Hoog, Knight & Stokes algorithm.
  
  The un-accelerated Fourier series approach is:
  
  .. math ::
  
      f(t,2M+1) = \frac{e^{\gamma t}}{T} \sum_{k=0}^{2M}{}^{'}
      \Re\left[\bar{f}\left( p_k \right)
      e^{i\pi t/T} \right],
  
  where the prime on the summation indicates the first term is halved.
  
  This simplistic approach requires so many function evaluations
  that it is not practical. Non-linear acceleration is
  accomplished via Pade-approximation and an analytic expression
  for the remainder of the continued fraction. See the original
  paper (reference 2 below) a detailed description of the
  numerical approach.
  
  **References**
  
  1. Davies, B. (2005). *Integral Transforms and their
     Applications*, Third Edition. Springer.
  2. de Hoog, F., J. Knight, A. Stokes (1982). An improved
     method for numerical inversion of Laplace transforms. *SIAM
     Journal of Scientific and Statistical Computing* 3:357-366,
     http://dx.doi.org/10.1137/0903022
docstr_tok: ["Calculate","time-domain","solution","for","de","Hoog",",","Knight","&","Stokes","algorithm",".","The","un-accelerated","Fourier","series","approach","is",":","..","math",":",":","f","(","t,2M+1",")","=","\\frac","{","e^","{","\\gamma","t","}","}","{","T","}","\\sum_","{","k=0","}","^","{","2M","}","{","}","^","{","'","}","\\Re\\left","[","\\bar","{","f","}","\\left","(","p_k","\\right",")","e^","{","i\\pi","t\/T","}","\\right","]",",","where","the","prime","on","the","summation","indicates","the","first","term","is","halved",".","This","simplistic","approach","requires","so","many","function","evaluations","that","it","is","not","practical",".","Non-linear","acceleration","is","accomplished","via","Pade-approximation","and","an","analytic","expression","for","the","remainder","of","the","continued","fraction",".","See","the","original","paper","(","reference","2","below",")","a","detailed","description","of","the","numerical","approach",".","*","*","References","*","*","1",".","Davies",",","B",".","(","2005",")",".","*","Integral","Transforms","and","their","Applications","*",",","Third","Edition",".","Springer",".","2.","de","Hoog",",","F.",",","J.","Knight",",","A.","Stokes","(","1982",")",".","An","improved","method","for","numerical","inversion","of","Laplace","transforms",".","*","SIAM","Journal","of","Scientific","and","Statistical","Computing","*","3:357-366",",","http",":","\/\/dx.doi.org\/10.1137\/0903022"]
code: def calc_time_domain_solution(self, fp, t, manual_prec=False):
      r"""Calculate time-domain solution for
      de Hoog, Knight & Stokes algorithm.
  
      The un-accelerated Fourier series approach is:
  
      .. math ::
  
          f(t,2M+1) = \frac{e^{\gamma t}}{T} \sum_{k=0}^{2M}{}^{'}
          \Re\left[\bar{f}\left( p_k \right)
          e^{i\pi t/T} \right],
  
      where the prime on the summation indicates the first term is halved.
  
      This simplistic approach requires so many function evaluations
      that it is not practical. Non-linear acceleration is
      accomplished via Pade-approximation and an analytic expression
      for the remainder of the continued fraction. See the original
      paper (reference 2 below) a detailed description of the
      numerical approach.
  
      **References**
  
      1. Davies, B. (2005). *Integral Transforms and their
         Applications*, Third Edition. Springer.
      2. de Hoog, F., J. Knight, A. Stokes (1982). An improved
         method for numerical inversion of Laplace transforms. *SIAM
         Journal of Scientific and Statistical Computing* 3:357-366,
         http://dx.doi.org/10.1137/0903022
  
      """
  
      M = self.degree
      np = self.np
      T = self.T
  
      self.t = self.ctx.convert(t)
  
      # would it be useful to try re-using
      # space between e&q and A&B?
      e = self.ctx.zeros(np, M + 1)
      q = self.ctx.matrix(2 * M, M)
      d = self.ctx.matrix(np, 1)
      A = self.ctx.zeros(np + 1, 1)
      B = self.ctx.ones(np + 1, 1)
  
      # initialize Q-D table
      e[:, 0] = 0.0 + 0j
      q[0, 0] = fp[1] / (fp[0] / 2)
      for i in range(1, 2 * M):
          q[i, 0] = fp[i + 1] / fp[i]
  
      # rhombus rule for filling triangular Q-D table (e & q)
      for r in range(1, M + 1):
          # start with e, column 1, 0:2*M-2
          mr = 2 * (M - r) + 1
          e[0:mr, r] = (
              q[1 : mr + 1, r - 1]
              - q[0:mr, r - 1]
              + e[1 : mr + 1, r - 1]
          )
          if not r == M:
              rq = r + 1
              mr = 2 * (M - rq) + 1 + 2
              for i in range(mr):
                  q[i, rq - 1] = (
                      q[i + 1, rq - 2] * e[i + 1, rq - 1] / e[i, rq - 1]
                  )
  
      # build up continued fraction coefficients (d)
      d[0] = fp[0] / 2
      for r in range(1, M + 1):
          d[2 * r - 1] = -q[0, r - 1]  # even terms
          d[2 * r] = -e[0, r]  # odd terms
  
      # seed A and B for recurrence
      A[0] = 0.0 + 0.0j
      A[1] = d[0]
      B[0:2] = 1.0 + 0.0j
  
      # base of the power series
      z = self.ctx.expjpi(self.t / T)  # i*pi is already in fcn
  
      # coefficients of Pade approximation (A & B)
      # using recurrence for all but last term
      for i in range(1, 2 * M):
          A[i + 1] = A[i] + d[i] * A[i - 1] * z
          B[i + 1] = B[i] + d[i] * B[i - 1] * z
  
      # "improved remainder" to continued fraction
      brem = (1 + (d[2 * M - 1] - d[2 * M]) * z) / 2
      # powm1(x,y) computes x^y - 1 more accurately near zero
      rem = brem * self.ctx.powm1(
          1 + d[2 * M] * z / brem, self.ctx.fraction(1, 2)
      )
  
      # last term of recurrence using new remainder
      A[np] = A[2 * M] + rem * A[2 * M - 1]
      B[np] = B[2 * M] + rem * B[2 * M - 1]
  
      # diagonal Pade approximation
      # F=A/B represents accelerated trapezoid rule
      result = (
          self.ctx.exp(self.gamma * self.t) / T * (A[np] / B[np]).real
      )
  
      # setting dps back to value when calc_laplace_parameter was called
      if not manual_prec:
          self.ctx.dps = self.dps_orig
  
      return result
  
label_desc: Use image node_id 2 for calling the deHoog obj's underlying member method code with example usage: obj.calc_time_domain_solution(fp, t, manual_prec) and returns: result