added trunc dir

Author: sveinlin

src/trunc/trunc_decay_FE.py

@@ -0,0 +1,25 @@
+"""
+Empirical estimation of the truncation error in a scheme.
+Examples on the Forward Euler scheme for the decay ODE u'=-au.
+"""
+import numpy as np
+import trunc_empir
+
+def decay_FE(dt, N):
+    dt = float(dt)
+    t = np.linspace(0, N*dt, N+1)
+    u_e = I*np.exp(-a*t)  # exact solution, I and a are global
+    u = u_e  # naming convention when writing up the scheme
+    R = np.zeros(N)
+
+    for n in range(0, N):
+        R[n] = (u[n+1] - u[n])/dt + a*u[n]
+
+    # Theoretical expression for the trunction error
+    R_a = 0.5*I*(-a)**2*np.exp(-a*t)*dt
+
+    return R, t[:-1], R_a[:-1]
+
+if __name__ == '__main__':
+    I = 1; a = 2  # global variables needed in decay_FE
+    trunc_empir.estimate(decay_FE, T=2.5, N_0=6, m=4, makeplot=True)

src/trunc/trunc_empir.py

@@ -0,0 +1,106 @@
+"""
+Empirical estimation of the truncation error in a scheme
+for a problem with one independent variable.
+"""
+import numpy as np
+import scitools.std as plt
+
+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], log='y')
+            legends_R.append('N=%d' % N[i])
+            plt.hold('on')
+
+            plt.figure(2)
+            plt.plot(t_coarse, R_a[i] - R[i], log='y')
+            plt.hold('on')
+            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:',
+    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')
+
+def convergence_rates(h, E):
+    """
+    Given a sequence of discretization parameters in the list h,
+    and corresponding errors in the list E,
+    compute the convergence rate of two successive (h[i], E[i])
+    and (h[i+1],E[i+1]) experiments, assuming the model E=C*h^r
+    (for small enough h).
+    """
+    from math import log
+    r = [log(E[i]/E[i-1])/log(h[i]/h[i-1])
+         for i in range(1, len(h))]
+    return r

src/trunc/truncation_errors.py

@@ -0,0 +1,102 @@
+import sympy as sym
+
+class TaylorSeries:
+    """Class for symbolic Taylor series."""
+    def __init__(self, f, num_terms=4):
+        self.f = f
+        self.N = num_terms
+        # Introduce symbols for the derivatives
+        self.df = [f]
+        for i in range(1, self.N+1):
+            self.df.append(sym.Symbol('D%d%s' % (i, f.name)))
+
+    def __call__(self, h):
+        """Return the truncated Taylor series at x+h."""
+        terms = self.f
+        for i in range(1, self.N+1):
+            terms += sym.Rational(1, sym.factorial(i))*self.df[i]*h**i
+        return terms
+
+
+class DiffOp:
+    """Class for discrete difference operators."""
+    def __init__(self, f, independent_variable='x',
+                 num_terms_Taylor_series=4):
+        self.Taylor = TaylorSeries(f, num_terms_Taylor_series)
+        self.f = self.Taylor.f
+        self.h = sym.Symbol('d%s' % independent_variable)
+
+        # Finite difference operators
+        h, f, f_T = self.h, self.f, self.Taylor # short names
+        theta = sym.Symbol('theta')
+        self.diffops = {
+            'Dtp':  (f_T(h) - f)/h,
+            'Dtm':  (f - f_T(-h))/h,
+            'Dt':   (f_T(h/2) - f_T(-h/2))/h,
+            'D2t':  (f_T(h) - f_T(-h))/(2*h),
+            'DtDt': (f_T(h) - 2*f + f_T(-h))/h**2,
+            'barDt': (f_T((1-theta)*h) - f_T(-theta*h))/h,
+            }
+        self.diffops = {diffop: sym.simplify(self.diffops[diffop])
+                        for diffop in self.diffops}
+
+        self.diffops['weighted_arithmetic_mean'] = \
+             self._weighted_arithmetic_mean()
+        self.diffops['geometric_mean'] = self._geometric_mean()
+        self.diffops['harmonic_mean'] = self._harmonic_mean()
+
+    def _weighted_arithmetic_mean(self):
+        # The expansion is around n*h + theta*h
+        h, f, f_T = self.h, self.f, self.Taylor
+        theta = sym.Symbol('theta')
+        f_n = f_T(-h*theta)
+        f_np1 = f_T((1-theta)*h)
+        a_mean = theta*f_np1 + (1-theta)*f_n
+        return sym.expand(a_mean)
+
+    def _geometric_mean(self):
+        h, f, f_T = self.h, self.f, self.Taylor
+        f_nmhalf = f_T(-h/2)
+        f_nphalf = f_T(h/2)
+        g_mean = f_nmhalf*f_nphalf
+        return sym.expand(g_mean)
+
+    def _harmonic_mean(self):
+        h, f, f_T = self.h, self.f, self.Taylor
+        f_nmhalf = f_T(-h/2)
+        f_nphalf = f_T(h/2)
+        h_mean = 2/(1/f_nmhalf + 1/f_nphalf)
+        return sym.expand(h_mean)
+
+    def D(self, i):
+        """Return the symbol for the i-th derivative."""
+        return self.Taylor.df[i]
+
+    def __getitem__(self, operator_name):
+        return self.diffops.get(operator_name, None)
+
+    def operator_names(self):
+        """Return all names for the operators."""
+        return list(self.diffops.keys())
+
+def truncation_errors():
+    # Make a table
+    u, theta = sym.symbols('u theta')
+    diffop = DiffOp(u, independent_variable='t',
+                    num_terms_Taylor_series=5)
+    D1u = diffop.D(1)   # symbol for du/dt
+    D2u = diffop.D(2)   # symbol for d^2u/dt^2
+    print 'R Dt:', diffop['Dt'] - D1u
+    print 'R Dtm:', diffop['Dtm'] - D1u
+    print 'R Dtp:', diffop['Dtp'] - D1u
+    print 'R barDt:', diffop['barDt'] - D1u
+    print 'R DtDt:', diffop['DtDt'] - D2u
+    print 'R weighted arithmetic mean:', diffop['weighted_arithmetic_mean'] - u
+    print 'R arithmetic mean:', diffop['weighted_arithmetic_mean'].subs(theta, sym.Rational(1,2)) - u
+    print 'R geometric mean:', diffop['geometric_mean'] - u
+    dt = diffop.h
+    print 'R harmonic mean:', (diffop['harmonic_mean'] - u).\
+          series(dt, 0, 3).as_leading_term(dt)
+
+if __name__ == '__main__':
+    truncation_errors()