[frames] | no frames]

# Source Code for Module hedge.timestep.stability

``` 1  """Automatic size finding for stability regions."""
2
3  from __future__ import division
4
6
8  This program is free software: you can redistribute it and/or modify
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
<!--
expandto(location.href);
// -->

```

 Generated by Epydoc 3.0.1 on Sat Aug 29 14:33:46 2009 http://epydoc.sourceforge.net