Some changes with approx

diff --git a/const.py b/const.py
index 5e00e9e..40d1476 100644 (file)
--- a/const.py
+++ b/const.py
@@ -24,3 +24,15 @@ def timestep(val):
     return power(10, double(c))
     
 HT = timestep(HY)
+
+def test_exp_stable(hy=HY, hz=HZ, ht=HZ, c=C):
+    gy = (c**2 * ht**2) / (hy**2)
+    gz = (c**2 * ht**2) / (hz**2)
+    return (gy + gz < 1)
+
+def test_imp_stable(hy=HY, ht=HT, c=C):
+    gy = (c**2 * ht**2) / (hy**2)
+    return (gy <= 1)
+    
+    
+    

diff --git a/err.py b/err.py
new file mode 100755 (executable)
index 0000000..cfe49aa
--- /dev/null
+++ b/err.py
@@ -0,0 +1,47 @@
+#!/usr/bin/python
+#-*- coding: utf-8 -*-
+
+from optparse import OptionParser
+
+import numpy as np
+
+import const as c
+from grid import ExplicitSolver, norm
+from waves import Exm
+
+
+def err_main():
+    '''
+    Calculates error (norm) series
+    '''
+    
+    # parser initialization
+    parser = OptionParser()
+    parser.add_option('--hy', dest='hy', default='1e-7')
+    parser.add_option('--hz', dest='hz', default='1e-7')
+    parser.add_option('--ht', dest='ht', default='1e-16')
+    (opts, args) = parser.parse_args()
+
+    # constants
+    time = np.double('1e-14')
+    ly = np.double('2e-5')
+    lz = np.double('2e-5')
+    
+    ht = np.double(opts.ht)
+    hy = np.double(opts.hy)
+    hz = np.double(opts.hz)
+
+    solver = ExplicitSolver(int(ly / hy), int(lz / hz), time, ht, solver='inline')
+    uh = solver.solve()
+
+    u = uh.copy()
+    M = 550
+    for i in xrange(u.shape[0]):
+        for j in xrange(u.shape[1]):
+            u[i, j] = Exm(i * hy, j * hz, time, M)
+
+    print norm(u, uh)
+    
+
+if __name__ == '__main__':
+    err_main()

diff --git a/grid.py b/grid.py
index 583f16d..7a300b7 100644 (file)
--- a/grid.py
+++ b/grid.py
@@ -4,7 +4,7 @@ import numpy as np
 from scipy import weave
 from scipy.weave import blitz, inline, converters
 from const import timestep
-from sweep import sweep
+from sweep import native_sweep, sweep
 
 def norm(lhs, rhs):
     if not (isinstance(lhs, np.ndarray)
@@ -14,14 +14,14 @@ def norm(lhs, rhs):
     if lhs.shape != rhs.shape:
         raise ValueError("arrays must be the same shape")
 
-    return np.max(np.abc(lhs - rhs))
+    return np.max(np.abs(lhs - rhs))
 
 
 class ExplicitSolver(object):
     '''Represents equation solver with time stepper    
     '''
 
-    def __init__(self, I, J, time, ht=None, ly=2e-5, lz=2e-5, lambda_=2e-6, C=3e8, solver='blitz', debug=False):
+    def __init__(self, I, J, time, ht=None, ly=2e-5, lz=2e-5, lambda_=2e-6, C=3e8, solver='inline', debug=False):
 
         if solver == 'python':
             self._step = self._python_step
@@ -199,7 +199,7 @@ class ImplicitSolver(object):
                 self._init_bc(2, k)
                 for i in xrange(1, self.I - 1):
                     f = self._make_rhs(i, k)
-                    self.u[2, i, 1:-1] = sweep(a, b, c, f)
+                    self.u[2, i, 1:-1] = native_sweep(a, b, c, f)
                 self.u[0] = self.u[1].copy()
                 self.u[1] = self.u[2].copy()
             self._solver = True

diff --git a/gui.py b/gui.py
index cbcb49d..d9a9231 100644 (file)
--- a/gui.py
+++ b/gui.py
@@ -22,6 +22,13 @@ from grid import ExplicitSolver, ImplicitSolver, norm
 from waves import Exm
 import const as c
 
+
+def split_number(x):
+    return ('%g' % x).split('e')
+
+def normal(x):
+    return np.float64('1e' + split_number(x)[1])
+
 class Dummy(object):
     pass
 
@@ -46,34 +53,12 @@ class AppForm(QMainWindow):
 
         self.splitter = QSplitter()
         self.splitter.addWidget(self.param_frame)
-        right_frame = QWidget()
-        right_layout = QVBoxLayout(right_frame)
-        right_layout.addWidget(self.work_frame)
         
-        scheme_frame = QWidget()
-        scheme_layout = QGridLayout(scheme_frame)
-        self.a_group = self._make_analytic()
-        self.exp_group = self._make_explicit()
-        self.imp_group = self._make_implicit()
-        scheme_layout.addWidget(self.a_group, 0, 0, 1, 1)
-        scheme_layout.addWidget(self.exp_group, 0, 1, 1, 1)
-        scheme_layout.addWidget(self.imp_group, 0, 2, 1, 1)
-
-        self.button_scheme_run = QPushButton(u'Пуск')
-        self.connect(self.button_scheme_run, SIGNAL("clicked()"), self.on_scheme_run)
-        self.connect(self.button_scheme_run, SIGNAL("clicked()"), self.on_deviation_run)
-        self.button_scheme_run.setFocus(Qt.ActiveWindowFocusReason)
-        
-        right_layout.addWidget(scheme_frame)
-        right_layout.addWidget(self.button_scheme_run)
-        
-        self.splitter.addWidget(right_frame)
+        self.splitter.addWidget(self.work_frame)
         self.main_layout = QGridLayout(self.main_frame)
         # self.main_layout.addWidget(self.param_frame, 0, 0)
         # self.main_layout.addWidget(self.work_frame, 0, 1)
         self.main_layout.addWidget(self.splitter)
-        
-
         self.setCentralWidget(self.main_frame)
 
     def init_param_frame(self):
@@ -139,6 +124,27 @@ class AppForm(QMainWindow):
         i_layout.addRow(u"Шаг по t", self.edit_imp_step_t)
         
         return i_group
+
+    def _make_norm_frame(self):
+        self.norm_frame = QWidget()
+        self.n_fig = Figure((5.0, 4.0), dpi = self.dpi)
+        self.n_canvas = FigureCanvas(self.n_fig)
+        self.n_canvas.setParent(self.main_frame)
+        self.n_axes = self.n_fig.add_subplot(111)
+        n_layout = QGridLayout(self.norm_frame)
+        n_layout.addWidget(self.n_canvas, 0, 0, 1, 3)
+        n_toolbar = NavigationToolbar(self.n_canvas, self.norm_frame)
+        n_layout.addWidget(n_toolbar, 1, 0, 1, 3)
+
+        ht_button = QPushButton("test ht")
+        self.connect(ht_button, SIGNAL("clicked()"), self.on_test_ht)
+        n_layout.addWidget(ht_button, 2, 0)
+        hy_button = QPushButton("test hy")
+        self.connect(hy_button, SIGNAL("clicked()"), self.on_test_hy)
+        n_layout.addWidget(hy_button, 2, 1)
+        hz_button = QPushButton("test hz")
+        self.connect(hz_button, SIGNAL("clicked()"), self.on_test_hz)
+        n_layout.addWidget(hz_button, 2, 2)
     
     def init_work_frame(self):
         self.work_frame = QTabWidget()
@@ -156,16 +162,24 @@ class AppForm(QMainWindow):
         self.mpl_toolbar = NavigationToolbar(self.z_canvas, self.graph_frame)
         # self.mpl_toolbar.setOrientation(Qt.Vertical)
         graph_layout.addWidget(self.mpl_toolbar, 1, 0, 1, 1)
+        
+        scheme_frame = QWidget()
+        scheme_layout = QGridLayout(scheme_frame)
+        self.a_group = self._make_analytic()
+        self.exp_group = self._make_explicit()
+        self.imp_group = self._make_implicit()
+        scheme_layout.addWidget(self.a_group, 0, 0, 1, 1)
+        scheme_layout.addWidget(self.exp_group, 0, 1, 1, 1)
+        scheme_layout.addWidget(self.imp_group, 0, 2, 1, 1)
 
-        self.norm_frame = QWidget()
-        self.n_fig = Figure((5.0, 4.0), dpi = self.dpi)
-        self.n_canvas = FigureCanvas(self.n_fig)
-        self.n_canvas.setParent(self.main_frame)
-        self.n_axes = self.n_fig.add_subplot(111)
-        n_layout = QGridLayout(self.norm_frame)
-        n_layout.addWidget(self.n_canvas, 0, 0, 1, 1)
-        n_toolbar = NavigationToolbar(self.n_canvas, self.norm_frame)
-        n_layout.addWidget(n_toolbar, 1, 0, 1, 1)
+        self.button_scheme_run = QPushButton(u'Пуск')
+        self.connect(self.button_scheme_run, SIGNAL("clicked()"), self.on_scheme_run)
+        self.button_scheme_run.setFocus(Qt.ActiveWindowFocusReason)
+
+        graph_layout.addWidget(scheme_frame, 2, 0)
+        graph_layout.addWidget(self.button_scheme_run, 3, 0)
+        
+        self._make_norm_frame()
         
         self.about_frame = QWidget()
 
@@ -173,6 +187,92 @@ class AppForm(QMainWindow):
         self.work_frame.addTab(self.norm_frame, u"&Аппроксимация")
         # self.work_frame.addTab(self.about_frame, u"&О программе")
 
+    def on_test_ht(self):
+        ly_str = unicode(self.edit_LY.text())
+        lz_str = unicode(self.edit_LZ.text())
+        ny_str = unicode(self.edit_dim_y.text())
+        nz_str = unicode(self.edit_dim_z.text())
+        time_str = unicode(self.edit_T.text())
+
+        M = 600
+
+        ly = np.float64(ly_str)
+        lz = np.float64(lz_str)
+        ny = int(ny_str)
+        nz = int(nz_str)
+        hy = (ly / ny)
+        hz = (lz / nz)
+        time = np.float64(time_str)
+
+        ht_exp = min(hy, hz) / (np.sqrt(2) * c.C)
+        ht = normal(ht_exp)
+
+        x_vals = []
+        y_vals = []
+        for n in xrange(1, 5):
+            print ht
+            x_vals.append(ht)
+            solver = ExplicitSolver(ny, nz, time, ht=ht)
+            u = solver.solve()
+            for i in xrange(u.shape[0]):
+                for j in xrange(u.shape[1]):
+                    u[i, j] -= Exm(i * hy, j * hz, time, M)
+            y_vals.append(np.max(np.abs(u)))
+            ht /= 2
+                    
+        line, = self.n_axes.plot(x_vals, y_vals, 'o-')
+
+        self.n_axes.grid(True)
+        self.n_axes.set_xlabel('ht')
+        self.n_axes.set_ylabel('err')
+        self.n_canvas.draw()
+
+    def on_test_hy(self):
+        ly_str = unicode(self.edit_LY.text())
+        lz_str = unicode(self.edit_LZ.text())
+        ny_str = unicode(self.edit_dim_y.text())
+        nz_str = unicode(self.edit_dim_z.text())
+        time_str = unicode(self.edit_T.text())
+
+        M = 600
+
+        ly = np.float64(ly_str)
+        lz = np.float64(lz_str)
+        ny = int(ny_str)
+        nz = int(nz_str)
+        hy = (ly / ny)
+        hz = (lz / nz)
+        time = np.float64(time_str)
+
+        ht_exp = min(hy, hz) / (np.sqrt(2) * c.C)
+        ht = normal(ht_exp)
+        ht /= 2
+        
+        x_vals = []
+        y_vals = []
+        hy_step = hy
+        for n in xrange(1, 11):
+            print hy
+            x_vals.append(hy)
+            solver = ExplicitSolver(ny, nz, time, ht=ht)
+            u = solver.solve()
+            for i in xrange(u.shape[0]):
+                for j in xrange(u.shape[1]):
+                    u[i, j] -= Exm(i * hy, j * hz, time, M)
+            y_vals.append(np.max(np.abs(u)))
+            hy /= 2
+            ny *= 2
+                    
+        line, = self.n_axes.plot(x_vals, y_vals, 'o-')
+
+        self.n_axes.grid(True)
+        self.n_axes.set_xlabel('ht')
+        self.n_axes.set_ylabel('err')
+        self.n_canvas.draw()
+
+    def on_test_hz(self):
+        pass
+
     def on_grid_change(self):
         ly_str = unicode(self.edit_LY.text())
         lz_str = unicode(self.edit_LZ.text())
@@ -228,9 +328,6 @@ class AppForm(QMainWindow):
         ly, lz, c, lambda_, t = [np.float64(unicode(p.text()).strip()) \
                                  for p in [self.edit_LY, self.edit_LZ, self.edit_C, self.edit_LAMBDA, self.edit_T]]
         return ly, lz, c, lambda_, t
-
-    def plot_graphs(self):
-        self.data.z
     
     def on_scheme_run(self):
         is_analytic, is_explicit, is_implicit = [g.isChecked() \
@@ -285,42 +382,42 @@ class AppForm(QMainWindow):
         self.z_axes.legend(lines, legs)
         self.z_canvas.draw()
 
-    def on_deviation_run(self):
-        flags = (is_analytic, is_explicit, is_implicit) = [g.isChecked() \
-                                                           for g in [self.a_group, self.exp_group, self.imp_group]]
-
-        if not (is_analytic or is_explicit or is_implicit):
-            return
-        if len(filter(None, flags)) < 2:
-            return
-
-        self.d_axes.clear()
-        lines = []
-        legs = []
-        if (is_analytic and is_explicit):
-            self.data.a_exp_dev = np.abs(self.data.E_a - self.data.E_exp)
-            line, = self.d_axes.plot(self.data.z, self.data.a_exp_dev, color='blue')
-            lines.append(line)
-            legs.append('Analytical and explicit scheme')
-        if (is_analytic and is_implicit):
-            self.data.a_imp_dev = np.abs(self.data.E_a - self.data.E_imp)
-            line, = self.d_axes.plot(self.data.z, self.data.a_imp_dev, color='green')
-            lines.append(line)
-            legs.append('Analytical and implicit scheme')
-        if (is_explicit and is_implicit):
-            self.data.exp_imp = np.abs(self.data.E_exp - self.data.E_imp)
-            line, = self.d_axes.plot(self.data.z, self.data.exp_imp, color='black')
-            lines.append(line)
-            legs.append('Explicit and implicit schemes')
-        self.d_axes.grid(True)
-        self.d_axes.set_xlabel('z')
-        self.d_axes.set_ylabel('e')
-        emax = max(self.data.E_a)
-        zmax = self.data.z[-1]
-        self.d_axes.set_ylim(0, 2*emax)
-        self.d_axes.set_xlim(0, zmax)
-        self.d_axes.legend(lines, legs)
-        self.d_canvas.draw()
+#     def on_deviation_run(self):
+#         flags = (is_analytic, is_explicit, is_implicit) = [g.isChecked() \
+#                                                            for g in [self.a_group, self.exp_group, self.imp_group]]
+
+#         if not (is_analytic or is_explicit or is_implicit):
+#             return
+#         if len(filter(None, flags)) < 2:
+#             return
+
+#         self.d_axes.clear()
+#         lines = []
+#         legs = []
+#         if (is_analytic and is_explicit):
+#             self.data.a_exp_dev = np.abs(self.data.E_a - self.data.E_exp)
+#             line, = self.d_axes.plot(self.data.z, self.data.a_exp_dev, color='blue')
+#             lines.append(line)
+#             legs.append('Analytical and explicit scheme')
+#         if (is_analytic and is_implicit):
+#             self.data.a_imp_dev = np.abs(self.data.E_a - self.data.E_imp)
+#             line, = self.d_axes.plot(self.data.z, self.data.a_imp_dev, color='green')
+#             lines.append(line)
+#             legs.append('Analytical and implicit scheme')
+#         if (is_explicit and is_implicit):
+#             self.data.exp_imp = np.abs(self.data.E_exp - self.data.E_imp)
+#             line, = self.d_axes.plot(self.data.z, self.data.exp_imp, color='black')
+#             lines.append(line)
+#             legs.append('Explicit and implicit schemes')
+#         self.d_axes.grid(True)
+#         self.d_axes.set_xlabel('z')
+#         self.d_axes.set_ylabel('e')
+#         emax = max(self.data.E_a)
+#         zmax = self.data.z[-1]
+#         self.d_axes.set_ylim(0, 2*emax)
+#         self.d_axes.set_xlim(0, zmax)
+#         self.d_axes.legend(lines, legs)
+#         self.d_canvas.draw()
 
 
 class SolverWorker(QThread):

diff --git a/sweep.py b/sweep.py
index 3d8ad39..2451a8f 100644 (file)
--- a/sweep.py
+++ b/sweep.py
@@ -2,34 +2,8 @@
 #-*- coding: utf-8 -*-
 
 import numpy as np
+from scipy.weave import blitz, inline, converters
 
-# def sweep(A, d):
-#     '''
-#     Решает систему Ax = d, где A - трехдиагональная матрица вида [0 ... a b c ... 0]
-#     d - вектор-строка
-    
-#     '''
-#     n = A.shape[0]
-#     a = np.zeros(n)
-#     a[1:] = np.diag(A, -1)
-#     b = np.diag(A)
-#     c = np.zeros(n)
-#     c[:-1] = np.diag(A, 1)
-#     x = d.copy()
-    
-#     # Прямой ход
-#     c[0] /= b[0]
-#     x[0] /= b[0]
-#     for i in xrange(1, n):
-#         denom = 1 / (b[i] - c[i-1] * a[i])
-#         c[i] *= denom
-#         x[i] = (x[i] - x[i-1] * a[i]) * denom
-
-#     # Обратный ход
-#     for i in xrange(n - 2, -1, -1):
-#         x[i] -= c[i] * x[i+1]
-
-#     return x
 
 def sweep(_a, _b, _c, _d):
     n = len(_b)
@@ -54,6 +28,33 @@ def sweep(_a, _b, _c, _d):
 
     return x
 
+def native_sweep(_a, _b, _c, _d):
+    n = len(_b)
+    a = np.zeros(n)
+    a[1:] = _a
+    b = _b.copy()
+    c = np.zeros(n)
+    c[:-1] = _c
+    x = _d.copy()
+
+    code = """
+    c(0) /= b(0);
+    x(0) /= b(0);
+    for (int i=1; i < n; i++) {
+        double denom = 1 / (b(i) - c(i - 1) * a(i));
+        c(i) *= denom;
+        x(i) = (x(i) - x(i - 1) * a(i)) * denom;
+    }
+    for (int i=n-2; i >= 0; i--) {
+        x(i) -= c(i) * x(i+1);
+    }
+    """
+    inline(code,
+           ['n', 'a', 'b', 'c', 'x'],
+           type_converters=converters.blitz,
+           compiler='gcc')
+    return x
+
 def main():
     pass