original_docstring: 
      Compute the truncation error in a problem with one independent
      variable, using m meshes, and estimate the convergence
      rate of the truncation error.
  
      The user-supplied function truncation_error(dt, N) computes
      the truncation error on a uniform mesh with N intervals of
      length dt::
  
        R, t, R_a = truncation_error(dt, N)
  
      where R holds the truncation error at points in the array t,
      and R_a are the corresponding theoretical truncation error
      values (None if not available).
  
      The truncation_error function is run on a series of meshes
      with 2**i*N_0 intervals, i=0,1,...,m-1.
      The values of R and R_a are restricted to the coarsest mesh.
      and based on these data, the convergence rate of R (pointwise)
      and time-integrated R can be estimated empirically.
      
docstring: Compute the truncation error in a problem with one independent
  variable, using m meshes, and estimate the convergence
  rate of the truncation error.
  
  The user-supplied function truncation_error(dt, N) computes
  the truncation error on a uniform mesh with N intervals of
  length dt:.
  
  
  
  
  
  The truncation_error function is run on a series of meshes
  with 2**i*N_0 intervals, i=0,1,...,m-1.
  The values of R and R_a are restricted to the coarsest mesh.
  and based on these data, the convergence rate of R (pointwise)
  and time-integrated R can be estimated empirically.
docstring_tokens: Compute the truncation error in a problem with one independent variable using m meshes and estimate the convergence rate of the truncation error . The user - supplied function truncation_error ( dt N ) computes the truncation error on a uniform mesh with N intervals of length dt : . The truncation_error function is run on a series of meshes with 2 ** i * N_0 intervals i = 0 1 ... m - 1 . The values of R and R_a are restricted to the coarsest mesh . and based on these data the convergence rate of R ( pointwise ) and time - integrated R can be estimated empirically .
code: def estimate(truncation_error, T, N_0, m, makeplot=True):
      """
      Compute the truncation error in a problem with one independent
      variable, using m meshes, and estimate the convergence
      rate of the truncation error.
  
      The user-supplied function truncation_error(dt, N) computes
      the truncation error on a uniform mesh with N intervals of
      length dt::
  
        R, t, R_a = truncation_error(dt, N)
  
      where R holds the truncation error at points in the array t,
      and R_a are the corresponding theoretical truncation error
      values (None if not available).
  
      The truncation_error function is run on a series of meshes
      with 2**i*N_0 intervals, i=0,1,...,m-1.
      The values of R and R_a are restricted to the coarsest mesh.
      and based on these data, the convergence rate of R (pointwise)
      and time-integrated R can be estimated empirically.
      """
      N = [2**i*N_0 for i in range(m)]
  
      R_I = np.zeros(m) # time-integrated R values on various meshes
      R   = [None]*m    # time series of R restricted to coarsest mesh
      R_a = [None]*m    # time series of R_a restricted to coarsest mesh
      dt = np.zeros(m)
      legends_R = [];  legends_R_a = []  # all legends of curves
  
      for i in range(m):
          dt[i] = T/float(N[i])
          R[i], t, R_a[i] = truncation_error(dt[i], N[i])
  
          R_I[i] = np.sqrt(dt[i]*np.sum(R[i]**2))
  
          if i == 0:
              t_coarse = t           # the coarsest mesh
  
          stride = N[i]//N_0
          R[i] = R[i][::stride]      # restrict to coarsest mesh
          R_a[i] = R_a[i][::stride]
  
          if makeplot:
              plt.figure(1)
              plt.plot(t_coarse, R[i])
              plt.yscale("log")
              legends_R.append('N=%d' % N[i])
  
              plt.figure(2)
              plt.plot(t_coarse, R_a[i] - R[i])
              plt.yscale("log")
              legends_R_a.append('N=%d' % N[i])
  
      if makeplot:
          plt.figure(1)
          plt.xlabel('time')
          plt.ylabel('pointwise truncation error')
          plt.legend(legends_R)
          plt.savefig('R_series.png')
          plt.savefig('R_series.pdf')
          plt.figure(2)
          plt.xlabel('time')
          plt.ylabel('pointwise error in estimated truncation error')
          plt.legend(legends_R_a)
          plt.savefig('R_error.png')
          plt.savefig('R_error.pdf')
  
      # Convergence rates
      r_R_I = convergence_rates(dt, R_I)
      print('R integrated in time; r:', end=' ')
      print(' '.join(['%.1f' % r for r in r_R_I]))
      R = np.array(R)  # two-dim. numpy array
      r_R = [convergence_rates(dt, R[:,n])[-1]
             for n in range(len(t_coarse))]
  
      # Plot convergence rates
      if makeplot:
          plt.figure()
          plt.plot(t_coarse, r_R)
          plt.xlabel('time')
          plt.ylabel('r')
          plt.axis([t_coarse[0], t_coarse[-1], 0, 2.5])
          plt.title('Pointwise rate $r$ in truncation error $\sim\Delta t^r$')
          plt.savefig('R_rate_series.png')
          plt.savefig('R_rate_series.pdf')
code_tokens: def estimate ( truncation_error , T , N_0 , m , makeplot = True ) : N = [ 2 ** i * N_0 for i in range ( m ) ] R_I = np . zeros ( m ) R = [ None ] * m R_a = [ None ] * m dt = np . zeros ( m ) legends_R = [ ] ; legends_R_a = [ ] for i in range ( m ) : dt [ i ] = T / float ( N [ i ] ) R [ i ] , t , R_a [ i ] = truncation_error ( dt [ i ] , N [ i ] ) R_I [ i ] = np . sqrt ( dt [ i ] * np . sum ( R [ i ] ** 2 ) ) if i == 0 : t_coarse = t stride = N [ i ] // N_0 R [ i ] = R [ i ] [ : : stride ] R_a [ i ] = R_a [ i ] [ : : stride ] if makeplot : plt . figure ( 1 ) plt . plot ( t_coarse , R [ i ] ) plt . yscale ( "log" ) legends_R . append ( 'N=%d' % N [ i ] ) plt . figure ( 2 ) plt . plot ( t_coarse , R_a [ i ] - R [ i ] ) plt . yscale ( "log" ) legends_R_a . append ( 'N=%d' % N [ i ] ) if makeplot : plt . figure ( 1 ) plt . xlabel ( 'time' ) plt . ylabel ( 'pointwise truncation error' ) plt . legend ( legends_R ) plt . savefig ( 'R_series.png' ) plt . savefig ( 'R_series.pdf' ) plt . figure ( 2 ) plt . xlabel ( 'time' ) plt . ylabel ( 'pointwise error in estimated truncation error' ) plt . legend ( legends_R_a ) plt . savefig ( 'R_error.png' ) plt . savefig ( 'R_error.pdf' ) r_R_I = convergence_rates ( dt , R_I ) print ( 'R integrated in time; r:' , end = ' ' ) print ( ' ' . join ( [ '%.1f' % r for r in r_R_I ] ) ) R = np . array ( R ) r_R = [ convergence_rates ( dt , R [ : , n ] ) [ - 1 ] for n in range ( len ( t_coarse ) ) ] if makeplot : plt . figure ( ) plt . plot ( t_coarse , r_R ) plt . xlabel ( 'time' ) plt . ylabel ( 'r' ) plt . axis ( [ t_coarse [ 0 ] , t_coarse [ - 1 ] , 0 , 2.5 ] ) plt . title ( 'Pointwise rate $r$ in truncation error $\sim\Delta t^r$' ) plt . savefig ( 'R_rate_series.png' ) plt . savefig ( 'R_rate_series.pdf' )
short_docstring: Compute the truncation error in a problem with one independent
  variable, using m meshes, and estimate the convergence
  rate of the truncation error.
short_docstring_tokens: Compute the truncation error in a problem with one independent variable using m meshes and estimate the convergence rate of the truncation error .
comment: """
      Compute the truncation error in a problem with one independent
      variable, using m meshes, and estimate the convergence
      rate of the truncation error.
  
      The user-supplied function truncation_error(dt, N) computes
      the truncation error on a uniform mesh with N intervals of
      length dt::
  
        R, t, R_a = truncation_error(dt, N)
  
      where R holds the truncation error at points in the array t,
      and R_a are the corresponding theoretical truncation error
      values (None if not available).
  
      The truncation_error function is run on a series of meshes
      with 2**i*N_0 intervals, i=0,1,...,m-1.
      The values of R and R_a are restricted to the coarsest mesh.
      and based on these data, the convergence rate of R (pointwise)
      and time-integrated R can be estimated empirically.
      """$SEP$# time-integrated R values on various meshes$SEP$# time series of R restricted to coarsest mesh$SEP$# time series of R_a restricted to coarsest mesh$SEP$# all legends of curves$SEP$# the coarsest mesh$SEP$# restrict to coarsest mesh$SEP$# Convergence rates$SEP$# two-dim. numpy array$SEP$# Plot convergence rates