"""
**Examples**::
>>> import nodepy.linear_multistep_method as lm
>>> ab3=lm.Adams_Bashforth(3)
>>> ab3.order()
3
>>> bdf2=lm.backward_difference_formula(2)
>>> bdf2.order()
2
>>> bdf2.is_zero_stable()
True
>>> bdf7=lm.backward_difference_formula(7)
>>> bdf7.is_zero_stable()
False
>>> bdf3=lm.backward_difference_formula(3)
>>> bdf3.A_alpha_stability()
86
>>> ssp32=lm.elm_ssp2(3)
>>> ssp32.order()
2
>>> ssp32.ssp_coefficient()
1/2
>>> ssp32.plot_stability_region() #doctest: +ELLIPSIS
<Figure...
"""
from __future__ import print_function
from __future__ import absolute_import
import numpy as np
import copy
import sympy
import nodepy.snp as snp
import matplotlib.pyplot as plt
from sympy import symbols, latex
from nodepy.snp import printable
from nodepy.general_linear_method import GeneralLinearMethod
from six.moves import map
from six.moves import range
from sympy import Rational
try:
import sympy.combinatorial as combinatorial
except ImportError:
import sympy.functions.combinatorial as combinatorial
[docs]class LinearMultistepMethod(GeneralLinearMethod):
r"""
Implementation of linear multistep methods in the form:
$$
\\alpha_k y_{n+k} + \\alpha_{k-1} y_{n+k-1} + \\cdots + \\alpha_0 y_n
= h ( \\beta_k f_{n+k} + \\cdots + \\beta_0 f_n )
$$
Methods are automatically normalized so that `\alpha_k=1`.
Notes: Representation follows Hairer & Wanner p. 368, NOT Butcher.
**References**:
* :cite:`hairer1993` Chapter III
* :cite:`butcher2003`
"""
def __init__(self,alpha,beta,name='Linear multistep method',shortname='LMM',
description=''):
self.beta = beta /alpha[-1]
self.alpha = alpha/alpha[-1]
self.name = name
self.shortname = shortname
if description != '':
self.info = description
else:
if self.is_explicit():
exp_str = "Explicit"
else:
exp_str = "Implicit"
self.info = "%s %d-step method of order %d" % (exp_str, len(self), self.order())
def __num__(self):
"""
Returns a copy of the method but with floating-point coefficients.
This is useful whenever we need to operate numerically without
worrying about the representation of the method.
"""
numself = copy.deepcopy(self)
if self.alpha.dtype==object:
numself.alpha=np.array(self.alpha,dtype=np.float64)
numself.beta=np.array(self.beta,dtype=np.float64)
return numself
def __str__(self):
"""
Pretty-prints the coefficients
"""
from nodepy.utils import array2strings
from nodepy.runge_kutta_method import _get_column_widths
alpha = array2strings(self.alpha, printzeros=True)
beta = array2strings(self.beta, printzeros=True)
_, colmax = _get_column_widths([alpha,beta])
s = self.name + '\n' + self.info + '\n'
s += 'alpha = ['
for alpha_i in alpha:
s += alpha_i.ljust(colmax + 1)
s += ']\n'
s += 'beta = ['
for beta_i in beta:
s += beta_i.ljust(colmax + 1)
s += ']'
return s.rstrip()
[docs] def characteristic_polynomials(self):
r"""
Returns the characteristic polynomials (also known as generating
polynomials) of a linear multistep method. They are:
$$ \\rho(z) = \\sum_{j=0}^k \\alpha_k z^k $$
$$ \\sigma(z) = \\sum_{j=0}^k \\beta_k z^k $$
**Examples**::
>>> from nodepy import lm
>>> ab5 = lm.Adams_Bashforth(5)
>>> rho,sigma = ab5.characteristic_polynomials()
>>> print(rho)
5 4
1 x - 1 x
>>> print(sigma)
4 3 2
2.64 x - 3.853 x + 3.633 x - 1.769 x + 0.3486
**Reference**: :cite:`hairer1993` p. 370, eq. 2.4
"""
rho=np.poly1d(self.alpha[::-1])
sigma=np.poly1d(self.beta[::-1])
return rho, sigma
[docs] def order(self,tol=1.e-10):
r"""Return the order of the local truncation error of a linear multistep method.
**Examples**::
>>> from nodepy import lm
>>> am3=lm.Adams_Moulton(3)
>>> am3.order()
4
"""
p = 0
while True:
if self._satisfies_order_conditions(p+1,tol):
p = p + 1
else:
return p
def _satisfies_order_conditions(self,p,tol):
""" Return *True* if the linear multistep method satisfies
the conditions of order `p` (only)."""
ii=snp.arange(len(self.alpha))
return abs(sum(ii**p*self.alpha-p*self.beta*ii**(p-1)))<tol
[docs] def absolute_monotonicity_radius(self):
return self.ssp_coefficient()
@property
def p(self):
return self.order()
[docs] def latex(self):
r"""
Print a LaTeX representation of a linear multistep formula.
**Example**::
>>> from nodepy import lm
>>> print(lm.Adams_Bashforth(2).latex())
\begin{align} y_{n + 2} - y_{n + 1} = \frac{3}{2}h f(y_{n + 1}) - \frac{1}{2}h f(y_{n})\end{align}
"""
n = symbols('n')
k = len(self)
alpha_terms = []
beta_terms = []
for i in range(k+1):
subscript = latex(n+k-i)
if self.alpha[k-i] != 0:
alpha_terms.append(printable(self.alpha[k-i],return_one=False) +
' y_{'+subscript+'}')
if self.beta[k-i] != 0:
beta_terms.append(printable(self.beta[k-i],return_one=False) +
'h f(y_{'+subscript+'})')
lhs = ' + '.join(alpha_terms)
rhs = ' + '.join(beta_terms)
s = r'\begin{align}'+ ' = '.join([lhs,rhs]) + r'\end{align}'
s = s.replace('+ -','-')
return s
[docs] def ssp_coefficient(self):
r"""Return the SSP coefficient of the method.
The SSP coefficient is given by
$$ \\min_{0 \\le j < k} - \\alpha_k/ \\beta_k $$
if `\alpha_j<0` and `\beta_j>0` for all `j`, and is equal to
zero otherwise.
**Examples**::
>>> from nodepy import lm
>>> ssp32=lm.elm_ssp2(3)
>>> ssp32.ssp_coefficient()
1/2
>>> bdf2=lm.backward_difference_formula(2)
>>> bdf2.ssp_coefficient()
0
"""
if np.any(self.alpha[:-1]>0) or np.any(self.beta<0):
return 0
return min([-self.alpha[j]/self.beta[j]
for j in range(len(self.alpha)-1) if self.beta[j]!=0])
[docs] def plot_stability_region(self,N=100,bounds=None,color='r',filled=True, alpha=1.,
to_file=False,longtitle=False):
r"""
The region of absolute stability of a linear multistep method is
the set
$$
\\{ z \\in C \\mid \\rho(\\zeta) - z \\sigma(\\zeta)
\\text{ satisfies the root condition} \\}
$$
where `\rho(\zeta)` and `\sigma(\zeta)` are the characteristic
functions of the method.
Also plots the boundary locus, which is
given by the set of points z:
$$
\\{z \\mid z =
\\rho(\\exp(\\imath \\theta)) /
\\sigma(\\exp(\\imath\\theta)), 0 \\le \\theta \\le 2\\pi \\}
$$
Here `\rho` and `\sigma` are the characteristic polynomials
of the method.
**Reference**: :cite:`leveque2007` section 7.6.1
**Input**: (all optional)
- N -- Number of gridpoints to use in each direction
- bounds -- limits of plotting region
- color -- color to use for this plot
- filled -- if true, stability region is filled in (solid); otherwise it is outlined
"""
rho, sigma = self.__num__().characteristic_polynomials()
mag = lambda z: _root_condition(rho-z*sigma)
vmag = np.vectorize(mag)
z = self._boundary_locus()
if bounds is None:
# Use boundary locus to decide plot region
realmax, realmin = np.max(np.real(z)), np.min(np.real(z))
imagmax, imagmin = np.max(np.imag(z)), np.min(np.imag(z))
deltar, deltai = realmax-realmin, imagmax-imagmin
bounds = (realmin-deltar/5.,realmax+deltar/5.,
imagmin-deltai/5.,imagmax+deltai/5.)
y=np.linspace(bounds[2],bounds[3],N)
Y=np.tile(y[:,np.newaxis],(1,N))
x=np.linspace(bounds[0],bounds[1],N)
X=np.tile(x,(N,1))
Z=X+Y*1j
R=1.5-vmag(Z)
if filled:
plt.contourf(X,Y,R,[0,1],colors=color,alpha=alpha)
else:
plt.contour(X,Y,R,[0,1],colors=color,alpha=alpha)
fig = plt.gcf()
ax = fig.get_axes()
if longtitle:
plt.setp(ax,title='Absolute Stability Region for '+self.name)
else:
plt.setp(ax,title='Stability region')
plt.plot([0,0],[bounds[2],bounds[3]],'--k',linewidth=2)
plt.plot([bounds[0],bounds[1]],[0,0],'--k',linewidth=2)
plt.plot(np.real(z),np.imag(z),color='k',linewidth=3)
plt.axis(bounds)
plt.axis('image')
if to_file:
plt.savefig(to_file, transparent=True, bbox_inches='tight', pad_inches=0.3)
plt.draw()
return fig
[docs] def plot_boundary_locus(self,N=1000,figsize=None):
r"""Plot the boundary locus, which is
given by the set of points
$$
\\{z \\mid z =
\\rho(\\exp(\\imath\\theta)) /
\\sigma(\\exp(\\imath\\theta)), 0\\le \\theta \\le 2\\pi \\}
$$
where `\rho` and `\sigma` are the characteristic polynomials
of the method.
**Reference**: :cite:`leveque2007` section 7.6.1
"""
z = self._boundary_locus(N)
if figsize is None:
plt.figure()
else:
plt.figure(figsize=figsize)
plt.plot(np.real(z),np.imag(z),color='k',linewidth=3)
plt.axis('image')
bounds = plt.axis()
plt.plot([0,0],[bounds[2],bounds[3]],'--k',linewidth=2)
plt.plot([bounds[0],bounds[1]],[0,0],'--k',linewidth=2)
plt.title('Boundary locus for '+self.name)
plt.draw()
def _boundary_locus(self, N=1000):
r"""Compute the boundary locus, which is
given by the set of points
$$
\\{z \\mid z =
\\rho(\\exp(\\imath\\theta)) /
\\sigma(\\exp(\\imath\\theta)), 0 \\le \\theta \\le 2\\pi \\}
$$
where `\rho` and `\sigma` are the characteristic polynomials
of the method.
**Reference**: :cite:`leveque2007` section 7.6.1
"""
theta=np.linspace(0.,2*np.pi,N)
zeta = np.exp(theta*1j)
rho,sigma=self.__num__().characteristic_polynomials()
z = rho(zeta)/sigma(zeta)
return z
[docs] def A_alpha_stability(self, N=1000, tol=1.e-14):
r"""Angle of `A(\alpha)`-stability.
The result is given in degrees. The result is only accurate to
about 1 degree, so we round down.
**Examples**:
>>> from nodepy import lm
>>> bdf5 = lm.backward_difference_formula(5)
>>> bdf5.A_alpha_stability()
51
"""
from math import atan2, floor
z = self._boundary_locus(N)
rad = list(map(atan2,np.imag(z),np.real(z)))
rad = np.mod(np.array(rad),2*np.pi)
return min(int(floor(np.min(np.abs(np.where(np.real(z)<-tol,rad,1.e99)-np.pi))/np.pi*180)),90)
[docs] def is_explicit(self):
return self.beta[-1]==0
[docs] def is_zero_stable(self,tol=1.e-13):
r"""*True* if the method is zero-stable.
**Examples**::
>>> from nodepy import lm
>>> bdf5=lm.backward_difference_formula(5)
>>> bdf5.is_zero_stable()
True
"""
rho, sigma = self.characteristic_polynomials()
return _root_condition(rho,tol)
def __len__(self):
r"""Returns the number of steps used."""
return len(self.alpha)-1
[docs]class AdditiveLinearMultistepMethod(GeneralLinearMethod):
r"""
Method for solving equations of the form
$$ y'(t) = f(y) + g(y) $$
The method takes the form
$$
\\alpha_k y_{n+k} + \\alpha_{k-1} y_{n+k-1} + \\cdots + \\alpha_0 y_n
= h ( \\beta_k f_{n+k} + \\cdots + \\beta_0 f_n
+ \\gamma_k f_{n+k} + \\cdots + \\gamma_0 f_n )
$$
Methods are automatically normalized so that `\alpha_k=1`.
The usual reference for these is Ascher, Ruuth, and Whetton.
But we follow a different notation (as just described).
"""
def __init__(self, alpha, beta, gamma, name='Additive linear multistep method'):
self.beta = beta /alpha[-1]
self.gamma = gamma/alpha[-1]
self.alpha = alpha/alpha[-1]
self.name = name
self.method1 = LinearMultistepMethod(alpha, beta)
self.method2 = LinearMultistepMethod(alpha, gamma)
def __num__(self):
"""
Returns a copy of the method but with floating-point coefficients.
This is useful whenever we need to operate numerically without
worrying about the representation of the method.
"""
numself = copy.deepcopy(self)
if self.alpha.dtype == object:
numself.alpha = np.array(self.alpha, dtype=np.float64)
numself.beta = np.array(self.beta, dtype=np.float64)
numself.gamma = np.array(self.gamma, dtype=np.float64)
return numself
[docs] def order(self,tol=1.e-10):
r""" Return the order of the local truncation error of an additive
linear multistep method. The output is the minimum of the
order of the component methods.
"""
orders = []
for method in (self.method1,self.method2):
p = 0
while True:
if method._satisfies_order_conditions(p+1,tol):
p = p + 1
else:
orders.append(p)
break
return min(orders)
[docs] def plot_imex_stability_region(self,both_real=False,N=100,color='r',filled=True,
alpha=1.,fignum=None,bounds=[-10, 1, -5, 5]):
r"""
**Input**: (all optional)
- N -- Number of gridpoints to use in each direction
- bounds -- limits of plotting region
- color -- color to use for this plot
- filled -- if true, stability region is filled in (solid); otherwise it is outlined
"""
rho, sigma1 = self.method1.__num__().characteristic_polynomials()
rho, sigma2 = self.method2.__num__().characteristic_polynomials()
if both_real:
mag = lambda a, b: _max_root(rho - a*sigma1 - b*sigma2)
else:
mag = lambda a, b: _max_root(rho - a*sigma1 - 1j*b*sigma2)
vmag = np.vectorize(mag)
y = np.linspace(bounds[2],bounds[3],N)
Y = np.tile(y[:,np.newaxis],(1,N))
x = np.linspace(bounds[0],bounds[1],N)
X = np.tile(x,(N,1))
R = vmag(X,Y)
h = plt.figure(fignum)
if filled:
plt.contourf(X,Y,R,[0,1],colors=color,alpha=alpha)
else:
plt.contour(X,Y,R,[0,1],colors=color,alpha=alpha)
plt.contour(X,Y,R,np.linspace(0,1,10),colors='k')
plt.title('IMEX Stability Region for '+self.name)
plt.plot([0,0],[bounds[2],bounds[3]],'--k',linewidth=2)
plt.plot([bounds[0],bounds[1]],[0,0],'--k',linewidth=2)
plt.axis(bounds)
return h
[docs] def stiff_damping_factor(self,epsilon=1.e-7):
r"""
Return the magnitude of the largest root at `z = -\infty`.
This routine just computes a numerical approximation
to the true value (with absolute accuracy *epsilon*).
"""
rho, sigma1 = self.method1.__num__().characteristic_polynomials()
rho, sigma2 = self.method2.__num__().characteristic_polynomials()
mag = lambda a, b: _max_root(rho - a*sigma1 - 1j*b*sigma2)
f=[]
z=-1.
f.append(mag(z,0))
while True:
z = z*10.
f.append(mag(z,0))
if np.abs(f[-1]-f[-2]) < epsilon:
return f[-1]
if len(f)>100:
print(f)
raise Exception('Unable to compute stiff damping factor: slow convergence')
# ======================================================
def _max_root(p):
return max(np.abs(p.r))
def _root_condition(p,tol=1.e-13):
r"""*True* if the polynomial `p` has all roots inside
the unit circle and roots on the boundary of the unit circle
are simple.
**Examples**::
>>> from nodepy import lm
>>> p = np.poly1d((1,0.4,2,0.5))
>>> lm._root_condition(p)
False
"""
if max(np.abs(p.r))>(1+tol):
return False
mod_one_roots = [r for r in p.r if abs(abs(r)-1)<tol]
for i,r1 in enumerate(mod_one_roots):
for r2 in mod_one_roots[i+1:]:
if abs(r1-r2)<tol:
return False
return True
# ======================================================
# Families of multistep methods
# ======================================================
[docs]def Adams_Bashforth(k):
r"""
Construct the k-step, Adams-Bashforth method.
The methods are explicit and have order k.
They have the form:
$$ y_{n+1} = y_n + h \\sum_{j=0}^{k-1} \\beta_j f(y_{n-k+j+1}) $$
They are generated using equations (1.5) and (1.7) from
:cite:`hairer1993` III.1, along with the binomial expansion.
**Examples**::
>>> import nodepy.linear_multistep_method as lm
>>> ab3=lm.Adams_Bashforth(3)
>>> print(ab3)
3-step Adams-Bashforth
Explicit 3-step method of order 3
alpha = [ 0 0 -1 1 ]
beta = [ 5/12 -4/3 23/12 0 ]
>>> ab3.order()
3
**Reference**: :cite:`hairer1993`
"""
one = Rational(1,1)
alpha=snp.zeros(k+1)
beta=snp.zeros(k+1)
alpha[k]=one
alpha[k-1]=-one
gamma=snp.zeros(k)
gamma[0]=one
beta[k-1]=one
betaj=snp.zeros(k+1)
for j in range(1,k):
gamma[j]=one-sum(gamma[:j]/snp.arange(j+1,1,-1))
for i in range(0,j+1):
betaj[k-i-1]=(-one)**i*combinatorial.factorials.binomial(j,i)*gamma[j]
beta=beta+betaj
name=str(k)+'-step Adams-Bashforth'
return LinearMultistepMethod(alpha,beta,name=name,shortname='AB'+str(k))
[docs]def Nystrom(k):
r"""
Construct the `k`-step explicit Nystrom linear multistep method.
The methods are explicit and have order `k`.
They have the form:
$$y_{n+1} = y_{n-1} + h \\sum_{j=0}^{k-1} \\beta_j f(y_{n-k+j+1})$$
They are generated using equations (1.13) and (1.7) from
:cite:`hairer1993` III.1, along with the binomial expansion
and the relation in exercise 4 on p. 367.
Note that the term "Nystrom method" is also commonly used to refer
to a class of methods for second-order ODEs; those are NOT
the methods generated by this function.
**Examples**::
>>> import nodepy.linear_multistep_method as lm
>>> nys3=lm.Nystrom(6)
>>> nys3.order()
6
**Reference**: :cite:`hairer1993`
"""
one = Rational(1,1)
alpha = snp.zeros(k+1)
alpha[k] = one
alpha[k-2] = -one
beta = snp.zeros(k+1)
kappa = snp.zeros(k)
gamma = snp.zeros(k)
gamma[0] = one
kappa[0] = 2*one
beta[k-1] = 2*one
betaj = snp.zeros(k+1)
for j in range(1,k):
gamma[j] = one-sum(gamma[:j]/snp.arange(j+1,1,-1))
kappa[j] = 2 * gamma[j] - gamma[j-1]
for i in range(0,j+1):
betaj[k-i-1] = (-one)**i*combinatorial.factorials.binomial(j,i)*kappa[j]
beta = beta+betaj
name = str(k)+'-step Nystrom'
return LinearMultistepMethod(alpha,beta,name=name,shortname='Nys'+str(k))
[docs]def Adams_Moulton(k):
r"""
Construct the `k`-step, Adams-Moulton method.
The methods are implicit and have order `k+1`.
They have the form:
$$ y_{n+1} = y_n + h \\sum_{j=0}^{k} \\beta_j f(y_{n-k+j+1}) $$
They are generated using equation (1.9) and the equation in
Exercise 3 from Hairer & Wanner III.1, along with the binomial
expansion.
**Examples**::
>>> import nodepy.linear_multistep_method as lm
>>> am3=lm.Adams_Moulton(3)
>>> am3.order()
4
**Reference**: :cite:`hairer1993`
"""
alpha=snp.zeros(k+1)
beta=snp.zeros(k+1)
alpha[k]=1
alpha[k-1]=-1
gamma=snp.zeros(k+1)
gamma[0]=1
beta[k]=1
betaj=snp.zeros(k+1)
for j in range(1,k+1):
gamma[j]= -sum(gamma[:j]/snp.arange(j+1,1,-1))
for i in range(0,j+1):
betaj[k-i]=(-1)**i*combinatorial.factorials.binomial(j,i)*gamma[j]
beta=beta+betaj
name=str(k)+'-step Adams-Moulton'
return LinearMultistepMethod(alpha,beta,name=name,shortname='AM'+str(k))
[docs]def Milne_Simpson(k):
r"""
Construct the `k`-step, Milne-Simpson method.
The methods are implicit and (for `k \ge 3`) have order `k+1`.
They have the form:
$$ y_{n+1} = y_{n-1} + h \\sum_{j=0}^{k} \\beta_j f(y_{n-k+j+1}) $$
They are generated using equation (1.15), the equation in
Exercise 3, and the relation in exercise 4, all from Hairer & Wanner
III.1, along with the binomial expansion.
**Examples**::
>>> import nodepy.linear_multistep_method as lm
>>> ms3=lm.Milne_Simpson(3)
>>> ms3.order()
4
**Reference**: :cite:`hairer1993`
"""
alpha = snp.zeros(k+1)
beta = snp.zeros(k+1)
alpha[k] = 1
alpha[k-2] = -1
gamma = snp.zeros(k+1)
kappa = snp.zeros(k+1)
gamma[0] = 1
kappa[0] = 2
beta[k] = 2
betaj = snp.zeros(k+1)
for j in range(1,k+1):
gamma[j] = -sum(gamma[:j]/snp.arange(j+1,1,-1))
kappa[j] = 2 * gamma[j] - gamma[j-1]
for i in range(0,j+1):
betaj[k-i] = (-1)**i*combinatorial.factorials.binomial(j,i)*kappa[j]
beta = beta+betaj
name = str(k)+'-step Milne-Simpson'
return LinearMultistepMethod(alpha,beta,name=name,shortname='MS'+str(k))
[docs]def elm_ssp2(k):
r"""
Returns the optimal SSP `k`-step linear multistep method of order 2.
**Examples**::
>>> import nodepy.linear_multistep_method as lm
>>> lm10=lm.elm_ssp2(10)
>>> lm10.ssp_coefficient()
8/9
"""
alpha=snp.zeros(k+1)
beta=snp.zeros(k+1)
alpha[-1]=sympy.Rational(1,1)
alpha[0]=sympy.Rational(-1,(k-1)**2)
alpha[k-1]=sympy.Rational(-(k-1)**2+1,(k-1)**2)
beta[k-1]=sympy.Rational(k,k-1)
name='Optimal '+str(k)+'-step, 2nd order SSP method.'
return LinearMultistepMethod(alpha,beta,name=name)
[docs]def sand_cc(s):
r"""Construct Sand's circle-contractive method of order `p=2(s+1)`
that uses `2^s + 1` steps.
**Examples**::
>>> import nodepy.linear_multistep_method as lm
>>> cc4 = lm.sand_cc(4)
>>> cc4.order()
10
>>> cc4.ssp_coefficient()
1/8
**Reference**: :cite:`sand1986`
"""
one = sympy.Rational(1)
zero = sympy.Rational(0)
k = 2**s + 1
Jn = [k,k-1]
for i in range(1,s+1):
Jn.append(k-1-2**i)
alpha = snp.zeros(k+1)
beta = snp.zeros(k+1)
# This is inefficient
for j in Jn:
tau_product = one
tau_sum = zero
tau = [one/(j-i) for i in Jn if i!=j]
tau_product = np.prod(tau)
tau_sum = np.sum(tau)
beta[j] = tau_product**2
alpha[j] = 2*beta[j]*tau_sum
return LinearMultistepMethod(alpha,beta,'Sand circle-contractive')
[docs]def arw2(gam,c):
r"""Returns the second order IMEX additive multistep method based on the
parametrization in Section 3.2 of Ascher, Ruuth, & Whetton. The parameters
are *gam* and *c*. Known methods are obtained with the following values:
(1/2,0): CNAB
(1/2,1/8): MCNAB
(0,1): CNLF
(1,0): SBDF
**Examples**::
>>> from nodepy import lm
>>> import sympy
>>> CNLF = lm.arw2(0,sympy.Rational(1))
>>> CNLF.order()
2
>>> CNLF.method1.ssp_coefficient()
1
>>> CNLF.method2.ssp_coefficient()
0
>>> print(CNLF.stiff_damping_factor()) #doctest: +ELLIPSIS
0.999...
"""
half = sympy.Rational(1,2)
alpha = snp.array([gam-half,-2*gam,gam+half])
beta = snp.array([c/2,1-gam-c,gam+c/2]) # implicit part
gamma = snp.array([-gam,gam+1,0]) # explicit part
return AdditiveLinearMultistepMethod(alpha,beta,gamma,'ARW2('+str(gam)+','+str(c)+')')
[docs]def arw3(gam,theta,c):
r"""Returns the third order IMEX additive multistep method based on the
parametrization in Section 3.3 of Ascher, Ruuth, & Whetton. The parameters
are *gam*, *theta*, and *c*. Known methods are obtained with the following
values:
(1,0,0): SBDF3
Note that there is one sign error in the ARW paper; it is corrected here.
"""
half = sympy.Rational(1,2)
third = sympy.Rational(1,3)
alpha = snp.array([-half*gam**2+third/2, 3*half*gam**2+gam-1, -3*half*gam**2-2*gam+half-theta,
half*gam**2+gam+third+theta])
beta = snp.array([5*half/6*theta-c, (gam**2-gam)*half+3*c-4*theta*third, 1-gam**2-3*c+23*theta*third/4,
(gam**2+gam)*half+c])
gamma = snp.array([(gam**2+gam)*half+5*theta*third/4, -gam**2-2*gam-4*theta*third,
(gam**2+3*gam)*half+1+23*theta*third/4,0])
return AdditiveLinearMultistepMethod(alpha,beta,gamma,'ARW3('+str(gam)+','+str(theta)+','+str(c)+')')
[docs]def loadLMM(which='All'):
"""
Load a set of standard linear multistep methods for testing.
**Examples**::
>>> from nodepy import lm
>>> ebdf5 = lm.loadLMM('eBDF5')
>>> print(ebdf5)
eBDF 5
Explicit 5-step method of order 5
alpha = [ -12/137 75/137 -200/137 300/137 -300/137 1 ]
beta = [ 60/137 -300/137 600/137 -600/137 300/137 0 ]
>>> ebdf5.is_zero_stable()
True
"""
LM={}
# ================================================
alpha = snp.array([-12,75,-200,300,-300,137])/sympy.Rational(137,1)
beta = snp.array([60,-300,600,-600,300,0])/sympy.Rational(137,1)
LM['eBDF5'] = LinearMultistepMethod(alpha,beta,'eBDF 5')
# ================================================
theta = sympy.Rational(1,2)
alpha = snp.array([-1,1])
beta = snp.array([1-theta,theta])
gamma = snp.array([1,0])
LM['ET112'] = AdditiveLinearMultistepMethod(alpha,beta,gamma,'Euler-Theta')
# ================================================
if which=='All':
return LM
else:
return LM[which]
if __name__ == "__main__":
import doctest
doctest.testmod()