Source Code for Module hedge.timestep.stability

 1  """Automatic size finding for stability regions.""" 
 2   
 3  from __future__ import division 
 4   
 5  __copyright__ = "Copyright (C) 2009 Andreas Stock, Andreas Kloeckner" 
 6   
 7  __license__ = """ 
 8  This program is free software: you can redistribute it and/or modify 
 9  it under the terms of the GNU General Public License as published by 
10  the Free Software Foundation, either version 3 of the License, or 
11  (at your option) any later version. 
12   
13  This program is distributed in the hope that it will be useful, 
14  but WITHOUT ANY WARRANTY; without even the implied warranty of 
15  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
16  GNU General Public License for more details. 
17   
18  You should have received a copy of the GNU General Public License 
19  along with this program.  If not, see U{http://www.gnu.org/licenses/}. 
20  """ 
21   
22   
23   
24   
25  import numpy 
26  from pytools import memoize 
27   
28   
29   
30   
31  # bisection based method to find bounds of stability region on Imaginary axis only --- 
32 -def calculate_fudged_stability_region(stepper_class, *stepper_args): 
33      return calculate_stability_region(stepper_class, *stepper_args) \ 
34              * stepper_class.dt_fudge_factor 
35   
36   
37   
38   
39  @memoize 
40 -def calculate_stability_region(stepper_class, *stepper_args): 
41      def stepper_maker(): 
42          return stepper_class(*stepper_args) 
43   
44      prec = 1e-5 
45   
46      def is_stable(stepper, k): 
47          y = 1 
48          for i in range(20): 
49              if abs(y) > 2: 
50                  return False 
51              y = stepper(y, i, 1, lambda t, y: k*y) 
52          return True 
53   
54      def make_k(angle, mag): 
55          from cmath import exp 
56          return -prec+mag*exp(1j*angle) 
57   
58      def refine(stepper_maker, angle, stable, unstable): 
59          assert is_stable(stepper_maker(), make_k(angle, stable)) 
60          assert not is_stable(stepper_maker(), make_k(angle, unstable)) 
61          while abs(stable-unstable) > prec: 
62              mid = (stable+unstable)/2 
63              if is_stable(stepper_maker(), make_k(angle, mid)): 
64                  stable = mid 
65              else: 
66                  unstable = mid 
67          else: 
68              return stable 
69   
70      def find_stable_k(stepper_maker, angle): 
71          mag = 1 
72   
73          if is_stable(stepper_maker(), make_k(angle, mag)): 
74              mag *= 2 
75              while is_stable(stepper_maker(), make_k(angle, mag)): 
76                  mag *= 2 
77   
78                  if mag > 2**8: 
79                      return mag 
80              return refine(stepper_maker, angle, mag/2, mag) 
81          else: 
82              mag /= 2 
83              while not is_stable(stepper_maker(), make_k(angle, mag)): 
84                  mag /= 2 
85   
86                  if mag < prec: 
87                      return mag 
88              return refine(stepper_maker, angle, mag, mag*2) 
89   
90      points = [] 
91      from cmath import pi 
92      for angle in numpy.array([pi/2, 3/2*pi]): 
93          points.append(make_k(angle, find_stable_k(stepper_maker, angle))) 
94   
95      points = numpy.array(points) 
96   
97      return abs(points[0]) 
98