remove venv
This commit is contained in:
Tykayn
tykayn
parent
056387013d
commit
0680c7594e
13999 changed files with 0 additions and 2895688 deletions
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
Binary file not shown.
|
|
@ -1,211 +0,0 @@
|
|||
import pytest
|
||||
|
||||
import numpy as np
|
||||
from numpy.testing import assert_allclose
|
||||
|
||||
from scipy.integrate import quad_vec
|
||||
|
||||
from multiprocessing.dummy import Pool
|
||||
|
||||
|
||||
quadrature_params = pytest.mark.parametrize(
|
||||
'quadrature', [None, "gk15", "gk21", "trapezoid"])
|
||||
|
||||
|
||||
@quadrature_params
|
||||
def test_quad_vec_simple(quadrature):
|
||||
n = np.arange(10)
|
||||
def f(x):
|
||||
return x ** n
|
||||
for epsabs in [0.1, 1e-3, 1e-6]:
|
||||
if quadrature == 'trapezoid' and epsabs < 1e-4:
|
||||
# slow: skip
|
||||
continue
|
||||
|
||||
kwargs = dict(epsabs=epsabs, quadrature=quadrature)
|
||||
|
||||
exact = 2**(n+1)/(n + 1)
|
||||
|
||||
res, err = quad_vec(f, 0, 2, norm='max', **kwargs)
|
||||
assert_allclose(res, exact, rtol=0, atol=epsabs)
|
||||
|
||||
res, err = quad_vec(f, 0, 2, norm='2', **kwargs)
|
||||
assert np.linalg.norm(res - exact) < epsabs
|
||||
|
||||
res, err = quad_vec(f, 0, 2, norm='max', points=(0.5, 1.0), **kwargs)
|
||||
assert_allclose(res, exact, rtol=0, atol=epsabs)
|
||||
|
||||
res, err, *rest = quad_vec(f, 0, 2, norm='max',
|
||||
epsrel=1e-8,
|
||||
full_output=True,
|
||||
limit=10000,
|
||||
**kwargs)
|
||||
assert_allclose(res, exact, rtol=0, atol=epsabs)
|
||||
|
||||
|
||||
@quadrature_params
|
||||
def test_quad_vec_simple_inf(quadrature):
|
||||
def f(x):
|
||||
return 1 / (1 + np.float64(x) ** 2)
|
||||
|
||||
for epsabs in [0.1, 1e-3, 1e-6]:
|
||||
if quadrature == 'trapezoid' and epsabs < 1e-4:
|
||||
# slow: skip
|
||||
continue
|
||||
|
||||
kwargs = dict(norm='max', epsabs=epsabs, quadrature=quadrature)
|
||||
|
||||
res, err = quad_vec(f, 0, np.inf, **kwargs)
|
||||
assert_allclose(res, np.pi/2, rtol=0, atol=max(epsabs, err))
|
||||
|
||||
res, err = quad_vec(f, 0, -np.inf, **kwargs)
|
||||
assert_allclose(res, -np.pi/2, rtol=0, atol=max(epsabs, err))
|
||||
|
||||
res, err = quad_vec(f, -np.inf, 0, **kwargs)
|
||||
assert_allclose(res, np.pi/2, rtol=0, atol=max(epsabs, err))
|
||||
|
||||
res, err = quad_vec(f, np.inf, 0, **kwargs)
|
||||
assert_allclose(res, -np.pi/2, rtol=0, atol=max(epsabs, err))
|
||||
|
||||
res, err = quad_vec(f, -np.inf, np.inf, **kwargs)
|
||||
assert_allclose(res, np.pi, rtol=0, atol=max(epsabs, err))
|
||||
|
||||
res, err = quad_vec(f, np.inf, -np.inf, **kwargs)
|
||||
assert_allclose(res, -np.pi, rtol=0, atol=max(epsabs, err))
|
||||
|
||||
res, err = quad_vec(f, np.inf, np.inf, **kwargs)
|
||||
assert_allclose(res, 0, rtol=0, atol=max(epsabs, err))
|
||||
|
||||
res, err = quad_vec(f, -np.inf, -np.inf, **kwargs)
|
||||
assert_allclose(res, 0, rtol=0, atol=max(epsabs, err))
|
||||
|
||||
res, err = quad_vec(f, 0, np.inf, points=(1.0, 2.0), **kwargs)
|
||||
assert_allclose(res, np.pi/2, rtol=0, atol=max(epsabs, err))
|
||||
|
||||
def f(x):
|
||||
return np.sin(x + 2) / (1 + x ** 2)
|
||||
exact = np.pi / np.e * np.sin(2)
|
||||
epsabs = 1e-5
|
||||
|
||||
res, err, info = quad_vec(f, -np.inf, np.inf, limit=1000, norm='max', epsabs=epsabs,
|
||||
quadrature=quadrature, full_output=True)
|
||||
assert info.status == 1
|
||||
assert_allclose(res, exact, rtol=0, atol=max(epsabs, 1.5 * err))
|
||||
|
||||
|
||||
def test_quad_vec_args():
|
||||
def f(x, a):
|
||||
return x * (x + a) * np.arange(3)
|
||||
a = 2
|
||||
exact = np.array([0, 4/3, 8/3])
|
||||
|
||||
res, err = quad_vec(f, 0, 1, args=(a,))
|
||||
assert_allclose(res, exact, rtol=0, atol=1e-4)
|
||||
|
||||
|
||||
def _lorenzian(x):
|
||||
return 1 / (1 + x**2)
|
||||
|
||||
|
||||
@pytest.mark.fail_slow(10)
|
||||
def test_quad_vec_pool():
|
||||
f = _lorenzian
|
||||
res, err = quad_vec(f, -np.inf, np.inf, norm='max', epsabs=1e-4, workers=4)
|
||||
assert_allclose(res, np.pi, rtol=0, atol=1e-4)
|
||||
|
||||
with Pool(10) as pool:
|
||||
def f(x):
|
||||
return 1 / (1 + x ** 2)
|
||||
res, _ = quad_vec(f, -np.inf, np.inf, norm='max', epsabs=1e-4, workers=pool.map)
|
||||
assert_allclose(res, np.pi, rtol=0, atol=1e-4)
|
||||
|
||||
|
||||
def _func_with_args(x, a):
|
||||
return x * (x + a) * np.arange(3)
|
||||
|
||||
|
||||
@pytest.mark.fail_slow(10)
|
||||
@pytest.mark.parametrize('extra_args', [2, (2,)])
|
||||
@pytest.mark.parametrize('workers', [1, 10])
|
||||
def test_quad_vec_pool_args(extra_args, workers):
|
||||
f = _func_with_args
|
||||
exact = np.array([0, 4/3, 8/3])
|
||||
|
||||
res, err = quad_vec(f, 0, 1, args=extra_args, workers=workers)
|
||||
assert_allclose(res, exact, rtol=0, atol=1e-4)
|
||||
|
||||
with Pool(workers) as pool:
|
||||
res, err = quad_vec(f, 0, 1, args=extra_args, workers=pool.map)
|
||||
assert_allclose(res, exact, rtol=0, atol=1e-4)
|
||||
|
||||
|
||||
@quadrature_params
|
||||
def test_num_eval(quadrature):
|
||||
def f(x):
|
||||
count[0] += 1
|
||||
return x**5
|
||||
|
||||
count = [0]
|
||||
res = quad_vec(f, 0, 1, norm='max', full_output=True, quadrature=quadrature)
|
||||
assert res[2].neval == count[0]
|
||||
|
||||
|
||||
def test_info():
|
||||
def f(x):
|
||||
return np.ones((3, 2, 1))
|
||||
|
||||
res, err, info = quad_vec(f, 0, 1, norm='max', full_output=True)
|
||||
|
||||
assert info.success is True
|
||||
assert info.status == 0
|
||||
assert info.message == 'Target precision reached.'
|
||||
assert info.neval > 0
|
||||
assert info.intervals.shape[1] == 2
|
||||
assert info.integrals.shape == (info.intervals.shape[0], 3, 2, 1)
|
||||
assert info.errors.shape == (info.intervals.shape[0],)
|
||||
|
||||
|
||||
def test_nan_inf():
|
||||
def f_nan(x):
|
||||
return np.nan
|
||||
|
||||
def f_inf(x):
|
||||
return np.inf if x < 0.1 else 1/x
|
||||
|
||||
res, err, info = quad_vec(f_nan, 0, 1, full_output=True)
|
||||
assert info.status == 3
|
||||
|
||||
res, err, info = quad_vec(f_inf, 0, 1, full_output=True)
|
||||
assert info.status == 3
|
||||
|
||||
|
||||
@pytest.mark.parametrize('a,b', [(0, 1), (0, np.inf), (np.inf, 0),
|
||||
(-np.inf, np.inf), (np.inf, -np.inf)])
|
||||
def test_points(a, b):
|
||||
# Check that initial interval splitting is done according to
|
||||
# `points`, by checking that consecutive sets of 15 point (for
|
||||
# gk15) function evaluations lie between `points`
|
||||
|
||||
points = (0, 0.25, 0.5, 0.75, 1.0)
|
||||
points += tuple(-x for x in points)
|
||||
|
||||
quadrature_points = 15
|
||||
interval_sets = []
|
||||
count = 0
|
||||
|
||||
def f(x):
|
||||
nonlocal count
|
||||
|
||||
if count % quadrature_points == 0:
|
||||
interval_sets.append(set())
|
||||
|
||||
count += 1
|
||||
interval_sets[-1].add(float(x))
|
||||
return 0.0
|
||||
|
||||
quad_vec(f, a, b, points=points, quadrature='gk15', limit=0)
|
||||
|
||||
# Check that all point sets lie in a single `points` interval
|
||||
for p in interval_sets:
|
||||
j = np.searchsorted(sorted(points), tuple(p))
|
||||
assert np.all(j == j[0])
|
||||
|
|
@ -1,305 +0,0 @@
|
|||
import itertools
|
||||
import pytest
|
||||
import numpy as np
|
||||
from numpy.testing import assert_allclose
|
||||
from scipy.integrate import ode
|
||||
|
||||
|
||||
def _band_count(a):
|
||||
"""Returns ml and mu, the lower and upper band sizes of a."""
|
||||
nrows, ncols = a.shape
|
||||
ml = 0
|
||||
for k in range(-nrows+1, 0):
|
||||
if np.diag(a, k).any():
|
||||
ml = -k
|
||||
break
|
||||
mu = 0
|
||||
for k in range(nrows-1, 0, -1):
|
||||
if np.diag(a, k).any():
|
||||
mu = k
|
||||
break
|
||||
return ml, mu
|
||||
|
||||
|
||||
def _linear_func(t, y, a):
|
||||
"""Linear system dy/dt = a * y"""
|
||||
return a.dot(y)
|
||||
|
||||
|
||||
def _linear_jac(t, y, a):
|
||||
"""Jacobian of a * y is a."""
|
||||
return a
|
||||
|
||||
|
||||
def _linear_banded_jac(t, y, a):
|
||||
"""Banded Jacobian."""
|
||||
ml, mu = _band_count(a)
|
||||
bjac = [np.r_[[0] * k, np.diag(a, k)] for k in range(mu, 0, -1)]
|
||||
bjac.append(np.diag(a))
|
||||
for k in range(-1, -ml-1, -1):
|
||||
bjac.append(np.r_[np.diag(a, k), [0] * (-k)])
|
||||
return bjac
|
||||
|
||||
|
||||
def _solve_linear_sys(a, y0, tend=1, dt=0.1,
|
||||
solver=None, method='bdf', use_jac=True,
|
||||
with_jacobian=False, banded=False):
|
||||
"""Use scipy.integrate.ode to solve a linear system of ODEs.
|
||||
|
||||
a : square ndarray
|
||||
Matrix of the linear system to be solved.
|
||||
y0 : ndarray
|
||||
Initial condition
|
||||
tend : float
|
||||
Stop time.
|
||||
dt : float
|
||||
Step size of the output.
|
||||
solver : str
|
||||
If not None, this must be "vode", "lsoda" or "zvode".
|
||||
method : str
|
||||
Either "bdf" or "adams".
|
||||
use_jac : bool
|
||||
Determines if the jacobian function is passed to ode().
|
||||
with_jacobian : bool
|
||||
Passed to ode.set_integrator().
|
||||
banded : bool
|
||||
Determines whether a banded or full jacobian is used.
|
||||
If `banded` is True, `lband` and `uband` are determined by the
|
||||
values in `a`.
|
||||
"""
|
||||
if banded:
|
||||
lband, uband = _band_count(a)
|
||||
else:
|
||||
lband = None
|
||||
uband = None
|
||||
|
||||
if use_jac:
|
||||
if banded:
|
||||
r = ode(_linear_func, _linear_banded_jac)
|
||||
else:
|
||||
r = ode(_linear_func, _linear_jac)
|
||||
else:
|
||||
r = ode(_linear_func)
|
||||
|
||||
if solver is None:
|
||||
if np.iscomplexobj(a):
|
||||
solver = "zvode"
|
||||
else:
|
||||
solver = "vode"
|
||||
|
||||
r.set_integrator(solver,
|
||||
with_jacobian=with_jacobian,
|
||||
method=method,
|
||||
lband=lband, uband=uband,
|
||||
rtol=1e-9, atol=1e-10,
|
||||
)
|
||||
t0 = 0
|
||||
r.set_initial_value(y0, t0)
|
||||
r.set_f_params(a)
|
||||
r.set_jac_params(a)
|
||||
|
||||
t = [t0]
|
||||
y = [y0]
|
||||
while r.successful() and r.t < tend:
|
||||
r.integrate(r.t + dt)
|
||||
t.append(r.t)
|
||||
y.append(r.y)
|
||||
|
||||
t = np.array(t)
|
||||
y = np.array(y)
|
||||
return t, y
|
||||
|
||||
|
||||
def _analytical_solution(a, y0, t):
|
||||
"""
|
||||
Analytical solution to the linear differential equations dy/dt = a*y.
|
||||
|
||||
The solution is only valid if `a` is diagonalizable.
|
||||
|
||||
Returns a 2-D array with shape (len(t), len(y0)).
|
||||
"""
|
||||
lam, v = np.linalg.eig(a)
|
||||
c = np.linalg.solve(v, y0)
|
||||
e = c * np.exp(lam * t.reshape(-1, 1))
|
||||
sol = e.dot(v.T)
|
||||
return sol
|
||||
|
||||
|
||||
@pytest.mark.thread_unsafe
|
||||
def test_banded_ode_solvers():
|
||||
# Test the "lsoda", "vode" and "zvode" solvers of the `ode` class
|
||||
# with a system that has a banded Jacobian matrix.
|
||||
|
||||
# This test does not test the Jacobian evaluation (banded or not)
|
||||
# of "lsoda" due to the nonstiff nature of the equations.
|
||||
|
||||
t_exact = np.linspace(0, 1.0, 5)
|
||||
|
||||
# --- Real arrays for testing the "lsoda" and "vode" solvers ---
|
||||
|
||||
# lband = 2, uband = 1:
|
||||
a_real = np.array([[-0.6, 0.1, 0.0, 0.0, 0.0],
|
||||
[0.2, -0.5, 0.9, 0.0, 0.0],
|
||||
[0.1, 0.1, -0.4, 0.1, 0.0],
|
||||
[0.0, 0.3, -0.1, -0.9, -0.3],
|
||||
[0.0, 0.0, 0.1, 0.1, -0.7]])
|
||||
|
||||
# lband = 0, uband = 1:
|
||||
a_real_upper = np.triu(a_real)
|
||||
|
||||
# lband = 2, uband = 0:
|
||||
a_real_lower = np.tril(a_real)
|
||||
|
||||
# lband = 0, uband = 0:
|
||||
a_real_diag = np.triu(a_real_lower)
|
||||
|
||||
real_matrices = [a_real, a_real_upper, a_real_lower, a_real_diag]
|
||||
real_solutions = []
|
||||
|
||||
for a in real_matrices:
|
||||
y0 = np.arange(1, a.shape[0] + 1)
|
||||
y_exact = _analytical_solution(a, y0, t_exact)
|
||||
real_solutions.append((y0, t_exact, y_exact))
|
||||
|
||||
def check_real(idx, solver, meth, use_jac, with_jac, banded):
|
||||
a = real_matrices[idx]
|
||||
y0, t_exact, y_exact = real_solutions[idx]
|
||||
t, y = _solve_linear_sys(a, y0,
|
||||
tend=t_exact[-1],
|
||||
dt=t_exact[1] - t_exact[0],
|
||||
solver=solver,
|
||||
method=meth,
|
||||
use_jac=use_jac,
|
||||
with_jacobian=with_jac,
|
||||
banded=banded)
|
||||
assert_allclose(t, t_exact)
|
||||
assert_allclose(y, y_exact)
|
||||
|
||||
for idx in range(len(real_matrices)):
|
||||
p = [['vode', 'lsoda'], # solver
|
||||
['bdf', 'adams'], # method
|
||||
[False, True], # use_jac
|
||||
[False, True], # with_jacobian
|
||||
[False, True]] # banded
|
||||
for solver, meth, use_jac, with_jac, banded in itertools.product(*p):
|
||||
check_real(idx, solver, meth, use_jac, with_jac, banded)
|
||||
|
||||
# --- Complex arrays for testing the "zvode" solver ---
|
||||
|
||||
# complex, lband = 2, uband = 1:
|
||||
a_complex = a_real - 0.5j * a_real
|
||||
|
||||
# complex, lband = 0, uband = 0:
|
||||
a_complex_diag = np.diag(np.diag(a_complex))
|
||||
|
||||
complex_matrices = [a_complex, a_complex_diag]
|
||||
complex_solutions = []
|
||||
|
||||
for a in complex_matrices:
|
||||
y0 = np.arange(1, a.shape[0] + 1) + 1j
|
||||
y_exact = _analytical_solution(a, y0, t_exact)
|
||||
complex_solutions.append((y0, t_exact, y_exact))
|
||||
|
||||
def check_complex(idx, solver, meth, use_jac, with_jac, banded):
|
||||
a = complex_matrices[idx]
|
||||
y0, t_exact, y_exact = complex_solutions[idx]
|
||||
t, y = _solve_linear_sys(a, y0,
|
||||
tend=t_exact[-1],
|
||||
dt=t_exact[1] - t_exact[0],
|
||||
solver=solver,
|
||||
method=meth,
|
||||
use_jac=use_jac,
|
||||
with_jacobian=with_jac,
|
||||
banded=banded)
|
||||
assert_allclose(t, t_exact)
|
||||
assert_allclose(y, y_exact)
|
||||
|
||||
for idx in range(len(complex_matrices)):
|
||||
p = [['bdf', 'adams'], # method
|
||||
[False, True], # use_jac
|
||||
[False, True], # with_jacobian
|
||||
[False, True]] # banded
|
||||
for meth, use_jac, with_jac, banded in itertools.product(*p):
|
||||
check_complex(idx, "zvode", meth, use_jac, with_jac, banded)
|
||||
|
||||
# lsoda requires a stiffer problem to switch to stiff solver
|
||||
# Use the Robertson equation with surrounding trivial equations to make banded
|
||||
|
||||
def stiff_f(t, y):
|
||||
return np.array([
|
||||
y[0],
|
||||
-0.04 * y[1] + 1e4 * y[2] * y[3],
|
||||
0.04 * y[1] - 1e4 * y[2] * y[3] - 3e7 * y[2]**2,
|
||||
3e7 * y[2]**2,
|
||||
y[4]
|
||||
])
|
||||
|
||||
def stiff_jac(t, y):
|
||||
return np.array([
|
||||
[1, 0, 0, 0, 0],
|
||||
[0, -0.04, 1e4*y[3], 1e4*y[2], 0],
|
||||
[0, 0.04, -1e4 * y[3] - 3e7 * 2 * y[2], -1e4*y[2], 0],
|
||||
[0, 0, 3e7*2*y[2], 0, 0],
|
||||
[0, 0, 0, 0, 1]
|
||||
])
|
||||
|
||||
def banded_stiff_jac(t, y):
|
||||
return np.array([
|
||||
[0, 0, 0, 1e4*y[2], 0],
|
||||
[0, 0, 1e4*y[3], -1e4*y[2], 0],
|
||||
[1, -0.04, -1e4*y[3]-3e7*2*y[2], 0, 1],
|
||||
[0, 0.04, 3e7*2*y[2], 0, 0]
|
||||
])
|
||||
|
||||
@pytest.mark.thread_unsafe
|
||||
def test_banded_lsoda():
|
||||
# expected solution is given by problem with full jacobian
|
||||
tfull, yfull = _solve_robertson_lsoda(use_jac=True, banded=False)
|
||||
|
||||
for use_jac in [True, False]:
|
||||
t, y = _solve_robertson_lsoda(use_jac, True)
|
||||
assert_allclose(t, tfull)
|
||||
assert_allclose(y, yfull)
|
||||
|
||||
def _solve_robertson_lsoda(use_jac, banded):
|
||||
|
||||
if use_jac:
|
||||
if banded:
|
||||
jac = banded_stiff_jac
|
||||
else:
|
||||
jac = stiff_jac
|
||||
else:
|
||||
jac = None
|
||||
|
||||
if banded:
|
||||
lband = 1
|
||||
uband = 2
|
||||
else:
|
||||
lband = None
|
||||
uband = None
|
||||
|
||||
r = ode(stiff_f, jac)
|
||||
r.set_integrator('lsoda',
|
||||
lband=lband, uband=uband,
|
||||
rtol=1e-9, atol=1e-10,
|
||||
)
|
||||
t0 = 0
|
||||
dt = 1
|
||||
tend = 10
|
||||
y0 = np.array([1.0, 1.0, 0.0, 0.0, 1.0])
|
||||
r.set_initial_value(y0, t0)
|
||||
|
||||
t = [t0]
|
||||
y = [y0]
|
||||
while r.successful() and r.t < tend:
|
||||
r.integrate(r.t + dt)
|
||||
t.append(r.t)
|
||||
y.append(r.y)
|
||||
|
||||
# Ensure that the Jacobian was evaluated
|
||||
# iwork[12] has the number of Jacobian evaluations.
|
||||
assert r._integrator.iwork[12] > 0
|
||||
|
||||
t = np.array(t)
|
||||
y = np.array(y)
|
||||
return t, y
|
||||
|
|
@ -1,714 +0,0 @@
|
|||
import sys
|
||||
|
||||
try:
|
||||
from StringIO import StringIO
|
||||
except ImportError:
|
||||
from io import StringIO
|
||||
|
||||
import numpy as np
|
||||
from numpy.testing import (assert_, assert_array_equal, assert_allclose,
|
||||
assert_equal)
|
||||
from pytest import raises as assert_raises
|
||||
|
||||
from scipy.sparse import coo_matrix
|
||||
from scipy.special import erf
|
||||
from scipy.integrate._bvp import (modify_mesh, estimate_fun_jac,
|
||||
estimate_bc_jac, compute_jac_indices,
|
||||
construct_global_jac, solve_bvp)
|
||||
|
||||
import pytest
|
||||
|
||||
|
||||
def exp_fun(x, y):
|
||||
return np.vstack((y[1], y[0]))
|
||||
|
||||
|
||||
def exp_fun_jac(x, y):
|
||||
df_dy = np.empty((2, 2, x.shape[0]))
|
||||
df_dy[0, 0] = 0
|
||||
df_dy[0, 1] = 1
|
||||
df_dy[1, 0] = 1
|
||||
df_dy[1, 1] = 0
|
||||
return df_dy
|
||||
|
||||
|
||||
def exp_bc(ya, yb):
|
||||
return np.hstack((ya[0] - 1, yb[0]))
|
||||
|
||||
|
||||
def exp_bc_complex(ya, yb):
|
||||
return np.hstack((ya[0] - 1 - 1j, yb[0]))
|
||||
|
||||
|
||||
def exp_bc_jac(ya, yb):
|
||||
dbc_dya = np.array([
|
||||
[1, 0],
|
||||
[0, 0]
|
||||
])
|
||||
dbc_dyb = np.array([
|
||||
[0, 0],
|
||||
[1, 0]
|
||||
])
|
||||
return dbc_dya, dbc_dyb
|
||||
|
||||
|
||||
def exp_sol(x):
|
||||
return (np.exp(-x) - np.exp(x - 2)) / (1 - np.exp(-2))
|
||||
|
||||
|
||||
def sl_fun(x, y, p):
|
||||
return np.vstack((y[1], -p[0]**2 * y[0]))
|
||||
|
||||
|
||||
def sl_fun_jac(x, y, p):
|
||||
n, m = y.shape
|
||||
df_dy = np.empty((n, 2, m))
|
||||
df_dy[0, 0] = 0
|
||||
df_dy[0, 1] = 1
|
||||
df_dy[1, 0] = -p[0]**2
|
||||
df_dy[1, 1] = 0
|
||||
|
||||
df_dp = np.empty((n, 1, m))
|
||||
df_dp[0, 0] = 0
|
||||
df_dp[1, 0] = -2 * p[0] * y[0]
|
||||
|
||||
return df_dy, df_dp
|
||||
|
||||
|
||||
def sl_bc(ya, yb, p):
|
||||
return np.hstack((ya[0], yb[0], ya[1] - p[0]))
|
||||
|
||||
|
||||
def sl_bc_jac(ya, yb, p):
|
||||
dbc_dya = np.zeros((3, 2))
|
||||
dbc_dya[0, 0] = 1
|
||||
dbc_dya[2, 1] = 1
|
||||
|
||||
dbc_dyb = np.zeros((3, 2))
|
||||
dbc_dyb[1, 0] = 1
|
||||
|
||||
dbc_dp = np.zeros((3, 1))
|
||||
dbc_dp[2, 0] = -1
|
||||
|
||||
return dbc_dya, dbc_dyb, dbc_dp
|
||||
|
||||
|
||||
def sl_sol(x, p):
|
||||
return np.sin(p[0] * x)
|
||||
|
||||
|
||||
def emden_fun(x, y):
|
||||
return np.vstack((y[1], -y[0]**5))
|
||||
|
||||
|
||||
def emden_fun_jac(x, y):
|
||||
df_dy = np.empty((2, 2, x.shape[0]))
|
||||
df_dy[0, 0] = 0
|
||||
df_dy[0, 1] = 1
|
||||
df_dy[1, 0] = -5 * y[0]**4
|
||||
df_dy[1, 1] = 0
|
||||
return df_dy
|
||||
|
||||
|
||||
def emden_bc(ya, yb):
|
||||
return np.array([ya[1], yb[0] - (3/4)**0.5])
|
||||
|
||||
|
||||
def emden_bc_jac(ya, yb):
|
||||
dbc_dya = np.array([
|
||||
[0, 1],
|
||||
[0, 0]
|
||||
])
|
||||
dbc_dyb = np.array([
|
||||
[0, 0],
|
||||
[1, 0]
|
||||
])
|
||||
return dbc_dya, dbc_dyb
|
||||
|
||||
|
||||
def emden_sol(x):
|
||||
return (1 + x**2/3)**-0.5
|
||||
|
||||
|
||||
def undefined_fun(x, y):
|
||||
return np.zeros_like(y)
|
||||
|
||||
|
||||
def undefined_bc(ya, yb):
|
||||
return np.array([ya[0], yb[0] - 1])
|
||||
|
||||
|
||||
def big_fun(x, y):
|
||||
f = np.zeros_like(y)
|
||||
f[::2] = y[1::2]
|
||||
return f
|
||||
|
||||
|
||||
def big_bc(ya, yb):
|
||||
return np.hstack((ya[::2], yb[::2] - 1))
|
||||
|
||||
|
||||
def big_sol(x, n):
|
||||
y = np.ones((2 * n, x.size))
|
||||
y[::2] = x
|
||||
return x
|
||||
|
||||
|
||||
def big_fun_with_parameters(x, y, p):
|
||||
""" Big version of sl_fun, with two parameters.
|
||||
|
||||
The two differential equations represented by sl_fun are broadcast to the
|
||||
number of rows of y, rotating between the parameters p[0] and p[1].
|
||||
Here are the differential equations:
|
||||
|
||||
dy[0]/dt = y[1]
|
||||
dy[1]/dt = -p[0]**2 * y[0]
|
||||
dy[2]/dt = y[3]
|
||||
dy[3]/dt = -p[1]**2 * y[2]
|
||||
dy[4]/dt = y[5]
|
||||
dy[5]/dt = -p[0]**2 * y[4]
|
||||
dy[6]/dt = y[7]
|
||||
dy[7]/dt = -p[1]**2 * y[6]
|
||||
.
|
||||
.
|
||||
.
|
||||
|
||||
"""
|
||||
f = np.zeros_like(y)
|
||||
f[::2] = y[1::2]
|
||||
f[1::4] = -p[0]**2 * y[::4]
|
||||
f[3::4] = -p[1]**2 * y[2::4]
|
||||
return f
|
||||
|
||||
|
||||
def big_fun_with_parameters_jac(x, y, p):
|
||||
# big version of sl_fun_jac, with two parameters
|
||||
n, m = y.shape
|
||||
df_dy = np.zeros((n, n, m))
|
||||
df_dy[range(0, n, 2), range(1, n, 2)] = 1
|
||||
df_dy[range(1, n, 4), range(0, n, 4)] = -p[0]**2
|
||||
df_dy[range(3, n, 4), range(2, n, 4)] = -p[1]**2
|
||||
|
||||
df_dp = np.zeros((n, 2, m))
|
||||
df_dp[range(1, n, 4), 0] = -2 * p[0] * y[range(0, n, 4)]
|
||||
df_dp[range(3, n, 4), 1] = -2 * p[1] * y[range(2, n, 4)]
|
||||
|
||||
return df_dy, df_dp
|
||||
|
||||
|
||||
def big_bc_with_parameters(ya, yb, p):
|
||||
# big version of sl_bc, with two parameters
|
||||
return np.hstack((ya[::2], yb[::2], ya[1] - p[0], ya[3] - p[1]))
|
||||
|
||||
|
||||
def big_bc_with_parameters_jac(ya, yb, p):
|
||||
# big version of sl_bc_jac, with two parameters
|
||||
n = ya.shape[0]
|
||||
dbc_dya = np.zeros((n + 2, n))
|
||||
dbc_dyb = np.zeros((n + 2, n))
|
||||
|
||||
dbc_dya[range(n // 2), range(0, n, 2)] = 1
|
||||
dbc_dyb[range(n // 2, n), range(0, n, 2)] = 1
|
||||
|
||||
dbc_dp = np.zeros((n + 2, 2))
|
||||
dbc_dp[n, 0] = -1
|
||||
dbc_dya[n, 1] = 1
|
||||
dbc_dp[n + 1, 1] = -1
|
||||
dbc_dya[n + 1, 3] = 1
|
||||
|
||||
return dbc_dya, dbc_dyb, dbc_dp
|
||||
|
||||
|
||||
def big_sol_with_parameters(x, p):
|
||||
# big version of sl_sol, with two parameters
|
||||
return np.vstack((np.sin(p[0] * x), np.sin(p[1] * x)))
|
||||
|
||||
|
||||
def shock_fun(x, y):
|
||||
eps = 1e-3
|
||||
return np.vstack((
|
||||
y[1],
|
||||
-(x * y[1] + eps * np.pi**2 * np.cos(np.pi * x) +
|
||||
np.pi * x * np.sin(np.pi * x)) / eps
|
||||
))
|
||||
|
||||
|
||||
def shock_bc(ya, yb):
|
||||
return np.array([ya[0] + 2, yb[0]])
|
||||
|
||||
|
||||
def shock_sol(x):
|
||||
eps = 1e-3
|
||||
k = np.sqrt(2 * eps)
|
||||
return np.cos(np.pi * x) + erf(x / k) / erf(1 / k)
|
||||
|
||||
|
||||
def nonlin_bc_fun(x, y):
|
||||
# laplace eq.
|
||||
return np.stack([y[1], np.zeros_like(x)])
|
||||
|
||||
|
||||
def nonlin_bc_bc(ya, yb):
|
||||
phiA, phipA = ya
|
||||
phiC, phipC = yb
|
||||
|
||||
kappa, ioA, ioC, V, f = 1.64, 0.01, 1.0e-4, 0.5, 38.9
|
||||
|
||||
# Butler-Volmer Kinetics at Anode
|
||||
hA = 0.0-phiA-0.0
|
||||
iA = ioA * (np.exp(f*hA) - np.exp(-f*hA))
|
||||
res0 = iA + kappa * phipA
|
||||
|
||||
# Butler-Volmer Kinetics at Cathode
|
||||
hC = V - phiC - 1.0
|
||||
iC = ioC * (np.exp(f*hC) - np.exp(-f*hC))
|
||||
res1 = iC - kappa*phipC
|
||||
|
||||
return np.array([res0, res1])
|
||||
|
||||
|
||||
def nonlin_bc_sol(x):
|
||||
return -0.13426436116763119 - 1.1308709 * x
|
||||
|
||||
|
||||
def test_modify_mesh():
|
||||
x = np.array([0, 1, 3, 9], dtype=float)
|
||||
x_new = modify_mesh(x, np.array([0]), np.array([2]))
|
||||
assert_array_equal(x_new, np.array([0, 0.5, 1, 3, 5, 7, 9]))
|
||||
|
||||
x = np.array([-6, -3, 0, 3, 6], dtype=float)
|
||||
x_new = modify_mesh(x, np.array([1], dtype=int), np.array([0, 2, 3]))
|
||||
assert_array_equal(x_new, [-6, -5, -4, -3, -1.5, 0, 1, 2, 3, 4, 5, 6])
|
||||
|
||||
|
||||
def test_compute_fun_jac():
|
||||
x = np.linspace(0, 1, 5)
|
||||
y = np.empty((2, x.shape[0]))
|
||||
y[0] = 0.01
|
||||
y[1] = 0.02
|
||||
p = np.array([])
|
||||
df_dy, df_dp = estimate_fun_jac(lambda x, y, p: exp_fun(x, y), x, y, p)
|
||||
df_dy_an = exp_fun_jac(x, y)
|
||||
assert_allclose(df_dy, df_dy_an)
|
||||
assert_(df_dp is None)
|
||||
|
||||
x = np.linspace(0, np.pi, 5)
|
||||
y = np.empty((2, x.shape[0]))
|
||||
y[0] = np.sin(x)
|
||||
y[1] = np.cos(x)
|
||||
p = np.array([1.0])
|
||||
df_dy, df_dp = estimate_fun_jac(sl_fun, x, y, p)
|
||||
df_dy_an, df_dp_an = sl_fun_jac(x, y, p)
|
||||
assert_allclose(df_dy, df_dy_an)
|
||||
assert_allclose(df_dp, df_dp_an)
|
||||
|
||||
x = np.linspace(0, 1, 10)
|
||||
y = np.empty((2, x.shape[0]))
|
||||
y[0] = (3/4)**0.5
|
||||
y[1] = 1e-4
|
||||
p = np.array([])
|
||||
df_dy, df_dp = estimate_fun_jac(lambda x, y, p: emden_fun(x, y), x, y, p)
|
||||
df_dy_an = emden_fun_jac(x, y)
|
||||
assert_allclose(df_dy, df_dy_an)
|
||||
assert_(df_dp is None)
|
||||
|
||||
|
||||
def test_compute_bc_jac():
|
||||
ya = np.array([-1.0, 2])
|
||||
yb = np.array([0.5, 3])
|
||||
p = np.array([])
|
||||
dbc_dya, dbc_dyb, dbc_dp = estimate_bc_jac(
|
||||
lambda ya, yb, p: exp_bc(ya, yb), ya, yb, p)
|
||||
dbc_dya_an, dbc_dyb_an = exp_bc_jac(ya, yb)
|
||||
assert_allclose(dbc_dya, dbc_dya_an)
|
||||
assert_allclose(dbc_dyb, dbc_dyb_an)
|
||||
assert_(dbc_dp is None)
|
||||
|
||||
ya = np.array([0.0, 1])
|
||||
yb = np.array([0.0, -1])
|
||||
p = np.array([0.5])
|
||||
dbc_dya, dbc_dyb, dbc_dp = estimate_bc_jac(sl_bc, ya, yb, p)
|
||||
dbc_dya_an, dbc_dyb_an, dbc_dp_an = sl_bc_jac(ya, yb, p)
|
||||
assert_allclose(dbc_dya, dbc_dya_an)
|
||||
assert_allclose(dbc_dyb, dbc_dyb_an)
|
||||
assert_allclose(dbc_dp, dbc_dp_an)
|
||||
|
||||
ya = np.array([0.5, 100])
|
||||
yb = np.array([-1000, 10.5])
|
||||
p = np.array([])
|
||||
dbc_dya, dbc_dyb, dbc_dp = estimate_bc_jac(
|
||||
lambda ya, yb, p: emden_bc(ya, yb), ya, yb, p)
|
||||
dbc_dya_an, dbc_dyb_an = emden_bc_jac(ya, yb)
|
||||
assert_allclose(dbc_dya, dbc_dya_an)
|
||||
assert_allclose(dbc_dyb, dbc_dyb_an)
|
||||
assert_(dbc_dp is None)
|
||||
|
||||
|
||||
def test_compute_jac_indices():
|
||||
n = 2
|
||||
m = 4
|
||||
k = 2
|
||||
i, j = compute_jac_indices(n, m, k)
|
||||
s = coo_matrix((np.ones_like(i), (i, j))).toarray()
|
||||
s_true = np.array([
|
||||
[1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
|
||||
[1, 1, 1, 1, 0, 0, 0, 0, 1, 1],
|
||||
[0, 0, 1, 1, 1, 1, 0, 0, 1, 1],
|
||||
[0, 0, 1, 1, 1, 1, 0, 0, 1, 1],
|
||||
[0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
|
||||
[0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
|
||||
[1, 1, 0, 0, 0, 0, 1, 1, 1, 1],
|
||||
[1, 1, 0, 0, 0, 0, 1, 1, 1, 1],
|
||||
[1, 1, 0, 0, 0, 0, 1, 1, 1, 1],
|
||||
[1, 1, 0, 0, 0, 0, 1, 1, 1, 1],
|
||||
])
|
||||
assert_array_equal(s, s_true)
|
||||
|
||||
|
||||
def test_compute_global_jac():
|
||||
n = 2
|
||||
m = 5
|
||||
k = 1
|
||||
i_jac, j_jac = compute_jac_indices(2, 5, 1)
|
||||
x = np.linspace(0, 1, 5)
|
||||
h = np.diff(x)
|
||||
y = np.vstack((np.sin(np.pi * x), np.pi * np.cos(np.pi * x)))
|
||||
p = np.array([3.0])
|
||||
|
||||
f = sl_fun(x, y, p)
|
||||
|
||||
x_middle = x[:-1] + 0.5 * h
|
||||
y_middle = 0.5 * (y[:, :-1] + y[:, 1:]) - h/8 * (f[:, 1:] - f[:, :-1])
|
||||
|
||||
df_dy, df_dp = sl_fun_jac(x, y, p)
|
||||
df_dy_middle, df_dp_middle = sl_fun_jac(x_middle, y_middle, p)
|
||||
dbc_dya, dbc_dyb, dbc_dp = sl_bc_jac(y[:, 0], y[:, -1], p)
|
||||
|
||||
J = construct_global_jac(n, m, k, i_jac, j_jac, h, df_dy, df_dy_middle,
|
||||
df_dp, df_dp_middle, dbc_dya, dbc_dyb, dbc_dp)
|
||||
J = J.toarray()
|
||||
|
||||
def J_block(h, p):
|
||||
return np.array([
|
||||
[h**2*p**2/12 - 1, -0.5*h, -h**2*p**2/12 + 1, -0.5*h],
|
||||
[0.5*h*p**2, h**2*p**2/12 - 1, 0.5*h*p**2, 1 - h**2*p**2/12]
|
||||
])
|
||||
|
||||
J_true = np.zeros((m * n + k, m * n + k))
|
||||
for i in range(m - 1):
|
||||
J_true[i * n: (i + 1) * n, i * n: (i + 2) * n] = J_block(h[i], p[0])
|
||||
|
||||
J_true[:(m - 1) * n:2, -1] = p * h**2/6 * (y[0, :-1] - y[0, 1:])
|
||||
J_true[1:(m - 1) * n:2, -1] = p * (h * (y[0, :-1] + y[0, 1:]) +
|
||||
h**2/6 * (y[1, :-1] - y[1, 1:]))
|
||||
|
||||
J_true[8, 0] = 1
|
||||
J_true[9, 8] = 1
|
||||
J_true[10, 1] = 1
|
||||
J_true[10, 10] = -1
|
||||
|
||||
assert_allclose(J, J_true, rtol=1e-10)
|
||||
|
||||
df_dy, df_dp = estimate_fun_jac(sl_fun, x, y, p)
|
||||
df_dy_middle, df_dp_middle = estimate_fun_jac(sl_fun, x_middle, y_middle, p)
|
||||
dbc_dya, dbc_dyb, dbc_dp = estimate_bc_jac(sl_bc, y[:, 0], y[:, -1], p)
|
||||
J = construct_global_jac(n, m, k, i_jac, j_jac, h, df_dy, df_dy_middle,
|
||||
df_dp, df_dp_middle, dbc_dya, dbc_dyb, dbc_dp)
|
||||
J = J.toarray()
|
||||
assert_allclose(J, J_true, rtol=2e-8, atol=2e-8)
|
||||
|
||||
|
||||
def test_parameter_validation():
|
||||
x = [0, 1, 0.5]
|
||||
y = np.zeros((2, 3))
|
||||
assert_raises(ValueError, solve_bvp, exp_fun, exp_bc, x, y)
|
||||
|
||||
x = np.linspace(0, 1, 5)
|
||||
y = np.zeros((2, 4))
|
||||
assert_raises(ValueError, solve_bvp, exp_fun, exp_bc, x, y)
|
||||
|
||||
def fun(x, y, p):
|
||||
return exp_fun(x, y)
|
||||
def bc(ya, yb, p):
|
||||
return exp_bc(ya, yb)
|
||||
|
||||
y = np.zeros((2, x.shape[0]))
|
||||
assert_raises(ValueError, solve_bvp, fun, bc, x, y, p=[1])
|
||||
|
||||
def wrong_shape_fun(x, y):
|
||||
return np.zeros(3)
|
||||
|
||||
assert_raises(ValueError, solve_bvp, wrong_shape_fun, bc, x, y)
|
||||
|
||||
S = np.array([[0, 0]])
|
||||
assert_raises(ValueError, solve_bvp, exp_fun, exp_bc, x, y, S=S)
|
||||
|
||||
|
||||
def test_no_params():
|
||||
x = np.linspace(0, 1, 5)
|
||||
x_test = np.linspace(0, 1, 100)
|
||||
y = np.zeros((2, x.shape[0]))
|
||||
for fun_jac in [None, exp_fun_jac]:
|
||||
for bc_jac in [None, exp_bc_jac]:
|
||||
sol = solve_bvp(exp_fun, exp_bc, x, y, fun_jac=fun_jac,
|
||||
bc_jac=bc_jac)
|
||||
|
||||
assert_equal(sol.status, 0)
|
||||
assert_(sol.success)
|
||||
|
||||
assert_equal(sol.x.size, 5)
|
||||
|
||||
sol_test = sol.sol(x_test)
|
||||
|
||||
assert_allclose(sol_test[0], exp_sol(x_test), atol=1e-5)
|
||||
|
||||
f_test = exp_fun(x_test, sol_test)
|
||||
r = sol.sol(x_test, 1) - f_test
|
||||
rel_res = r / (1 + np.abs(f_test))
|
||||
norm_res = np.sum(rel_res**2, axis=0)**0.5
|
||||
assert_(np.all(norm_res < 1e-3))
|
||||
|
||||
assert_(np.all(sol.rms_residuals < 1e-3))
|
||||
assert_allclose(sol.sol(sol.x), sol.y, rtol=1e-10, atol=1e-10)
|
||||
assert_allclose(sol.sol(sol.x, 1), sol.yp, rtol=1e-10, atol=1e-10)
|
||||
|
||||
|
||||
def test_with_params():
|
||||
x = np.linspace(0, np.pi, 5)
|
||||
x_test = np.linspace(0, np.pi, 100)
|
||||
y = np.ones((2, x.shape[0]))
|
||||
|
||||
for fun_jac in [None, sl_fun_jac]:
|
||||
for bc_jac in [None, sl_bc_jac]:
|
||||
sol = solve_bvp(sl_fun, sl_bc, x, y, p=[0.5], fun_jac=fun_jac,
|
||||
bc_jac=bc_jac)
|
||||
|
||||
assert_equal(sol.status, 0)
|
||||
assert_(sol.success)
|
||||
|
||||
assert_(sol.x.size < 10)
|
||||
|
||||
assert_allclose(sol.p, [1], rtol=1e-4)
|
||||
|
||||
sol_test = sol.sol(x_test)
|
||||
|
||||
assert_allclose(sol_test[0], sl_sol(x_test, [1]),
|
||||
rtol=1e-4, atol=1e-4)
|
||||
|
||||
f_test = sl_fun(x_test, sol_test, [1])
|
||||
r = sol.sol(x_test, 1) - f_test
|
||||
rel_res = r / (1 + np.abs(f_test))
|
||||
norm_res = np.sum(rel_res ** 2, axis=0) ** 0.5
|
||||
assert_(np.all(norm_res < 1e-3))
|
||||
|
||||
assert_(np.all(sol.rms_residuals < 1e-3))
|
||||
assert_allclose(sol.sol(sol.x), sol.y, rtol=1e-10, atol=1e-10)
|
||||
assert_allclose(sol.sol(sol.x, 1), sol.yp, rtol=1e-10, atol=1e-10)
|
||||
|
||||
|
||||
def test_singular_term():
|
||||
x = np.linspace(0, 1, 10)
|
||||
x_test = np.linspace(0.05, 1, 100)
|
||||
y = np.empty((2, 10))
|
||||
y[0] = (3/4)**0.5
|
||||
y[1] = 1e-4
|
||||
S = np.array([[0, 0], [0, -2]])
|
||||
|
||||
for fun_jac in [None, emden_fun_jac]:
|
||||
for bc_jac in [None, emden_bc_jac]:
|
||||
sol = solve_bvp(emden_fun, emden_bc, x, y, S=S, fun_jac=fun_jac,
|
||||
bc_jac=bc_jac)
|
||||
|
||||
assert_equal(sol.status, 0)
|
||||
assert_(sol.success)
|
||||
|
||||
assert_equal(sol.x.size, 10)
|
||||
|
||||
sol_test = sol.sol(x_test)
|
||||
assert_allclose(sol_test[0], emden_sol(x_test), atol=1e-5)
|
||||
|
||||
f_test = emden_fun(x_test, sol_test) + S.dot(sol_test) / x_test
|
||||
r = sol.sol(x_test, 1) - f_test
|
||||
rel_res = r / (1 + np.abs(f_test))
|
||||
norm_res = np.sum(rel_res ** 2, axis=0) ** 0.5
|
||||
|
||||
assert_(np.all(norm_res < 1e-3))
|
||||
assert_allclose(sol.sol(sol.x), sol.y, rtol=1e-10, atol=1e-10)
|
||||
assert_allclose(sol.sol(sol.x, 1), sol.yp, rtol=1e-10, atol=1e-10)
|
||||
|
||||
|
||||
def test_complex():
|
||||
# The test is essentially the same as test_no_params, but boundary
|
||||
# conditions are turned into complex.
|
||||
x = np.linspace(0, 1, 5)
|
||||
x_test = np.linspace(0, 1, 100)
|
||||
y = np.zeros((2, x.shape[0]), dtype=complex)
|
||||
for fun_jac in [None, exp_fun_jac]:
|
||||
for bc_jac in [None, exp_bc_jac]:
|
||||
sol = solve_bvp(exp_fun, exp_bc_complex, x, y, fun_jac=fun_jac,
|
||||
bc_jac=bc_jac)
|
||||
|
||||
assert_equal(sol.status, 0)
|
||||
assert_(sol.success)
|
||||
|
||||
sol_test = sol.sol(x_test)
|
||||
|
||||
assert_allclose(sol_test[0].real, exp_sol(x_test), atol=1e-5)
|
||||
assert_allclose(sol_test[0].imag, exp_sol(x_test), atol=1e-5)
|
||||
|
||||
f_test = exp_fun(x_test, sol_test)
|
||||
r = sol.sol(x_test, 1) - f_test
|
||||
rel_res = r / (1 + np.abs(f_test))
|
||||
norm_res = np.sum(np.real(rel_res * np.conj(rel_res)),
|
||||
axis=0) ** 0.5
|
||||
assert_(np.all(norm_res < 1e-3))
|
||||
|
||||
assert_(np.all(sol.rms_residuals < 1e-3))
|
||||
assert_allclose(sol.sol(sol.x), sol.y, rtol=1e-10, atol=1e-10)
|
||||
assert_allclose(sol.sol(sol.x, 1), sol.yp, rtol=1e-10, atol=1e-10)
|
||||
|
||||
|
||||
def test_failures():
|
||||
x = np.linspace(0, 1, 2)
|
||||
y = np.zeros((2, x.size))
|
||||
res = solve_bvp(exp_fun, exp_bc, x, y, tol=1e-5, max_nodes=5)
|
||||
assert_equal(res.status, 1)
|
||||
assert_(not res.success)
|
||||
|
||||
x = np.linspace(0, 1, 5)
|
||||
y = np.zeros((2, x.size))
|
||||
res = solve_bvp(undefined_fun, undefined_bc, x, y)
|
||||
assert_equal(res.status, 2)
|
||||
assert_(not res.success)
|
||||
|
||||
|
||||
def test_big_problem():
|
||||
n = 30
|
||||
x = np.linspace(0, 1, 5)
|
||||
y = np.zeros((2 * n, x.size))
|
||||
sol = solve_bvp(big_fun, big_bc, x, y)
|
||||
|
||||
assert_equal(sol.status, 0)
|
||||
assert_(sol.success)
|
||||
|
||||
sol_test = sol.sol(x)
|
||||
|
||||
assert_allclose(sol_test[0], big_sol(x, n))
|
||||
|
||||
f_test = big_fun(x, sol_test)
|
||||
r = sol.sol(x, 1) - f_test
|
||||
rel_res = r / (1 + np.abs(f_test))
|
||||
norm_res = np.sum(np.real(rel_res * np.conj(rel_res)), axis=0) ** 0.5
|
||||
assert_(np.all(norm_res < 1e-3))
|
||||
|
||||
assert_(np.all(sol.rms_residuals < 1e-3))
|
||||
assert_allclose(sol.sol(sol.x), sol.y, rtol=1e-10, atol=1e-10)
|
||||
assert_allclose(sol.sol(sol.x, 1), sol.yp, rtol=1e-10, atol=1e-10)
|
||||
|
||||
|
||||
def test_big_problem_with_parameters():
|
||||
n = 30
|
||||
x = np.linspace(0, np.pi, 5)
|
||||
x_test = np.linspace(0, np.pi, 100)
|
||||
y = np.ones((2 * n, x.size))
|
||||
|
||||
for fun_jac in [None, big_fun_with_parameters_jac]:
|
||||
for bc_jac in [None, big_bc_with_parameters_jac]:
|
||||
sol = solve_bvp(big_fun_with_parameters, big_bc_with_parameters, x,
|
||||
y, p=[0.5, 0.5], fun_jac=fun_jac, bc_jac=bc_jac)
|
||||
|
||||
assert_equal(sol.status, 0)
|
||||
assert_(sol.success)
|
||||
|
||||
assert_allclose(sol.p, [1, 1], rtol=1e-4)
|
||||
|
||||
sol_test = sol.sol(x_test)
|
||||
|
||||
for isol in range(0, n, 4):
|
||||
assert_allclose(sol_test[isol],
|
||||
big_sol_with_parameters(x_test, [1, 1])[0],
|
||||
rtol=1e-4, atol=1e-4)
|
||||
assert_allclose(sol_test[isol + 2],
|
||||
big_sol_with_parameters(x_test, [1, 1])[1],
|
||||
rtol=1e-4, atol=1e-4)
|
||||
|
||||
f_test = big_fun_with_parameters(x_test, sol_test, [1, 1])
|
||||
r = sol.sol(x_test, 1) - f_test
|
||||
rel_res = r / (1 + np.abs(f_test))
|
||||
norm_res = np.sum(rel_res ** 2, axis=0) ** 0.5
|
||||
assert_(np.all(norm_res < 1e-3))
|
||||
|
||||
assert_(np.all(sol.rms_residuals < 1e-3))
|
||||
assert_allclose(sol.sol(sol.x), sol.y, rtol=1e-10, atol=1e-10)
|
||||
assert_allclose(sol.sol(sol.x, 1), sol.yp, rtol=1e-10, atol=1e-10)
|
||||
|
||||
|
||||
def test_shock_layer():
|
||||
x = np.linspace(-1, 1, 5)
|
||||
x_test = np.linspace(-1, 1, 100)
|
||||
y = np.zeros((2, x.size))
|
||||
sol = solve_bvp(shock_fun, shock_bc, x, y)
|
||||
|
||||
assert_equal(sol.status, 0)
|
||||
assert_(sol.success)
|
||||
|
||||
assert_(sol.x.size < 110)
|
||||
|
||||
sol_test = sol.sol(x_test)
|
||||
assert_allclose(sol_test[0], shock_sol(x_test), rtol=1e-5, atol=1e-5)
|
||||
|
||||
f_test = shock_fun(x_test, sol_test)
|
||||
r = sol.sol(x_test, 1) - f_test
|
||||
rel_res = r / (1 + np.abs(f_test))
|
||||
norm_res = np.sum(rel_res ** 2, axis=0) ** 0.5
|
||||
|
||||
assert_(np.all(norm_res < 1e-3))
|
||||
assert_allclose(sol.sol(sol.x), sol.y, rtol=1e-10, atol=1e-10)
|
||||
assert_allclose(sol.sol(sol.x, 1), sol.yp, rtol=1e-10, atol=1e-10)
|
||||
|
||||
|
||||
def test_nonlin_bc():
|
||||
x = np.linspace(0, 0.1, 5)
|
||||
x_test = x
|
||||
y = np.zeros([2, x.size])
|
||||
sol = solve_bvp(nonlin_bc_fun, nonlin_bc_bc, x, y)
|
||||
|
||||
assert_equal(sol.status, 0)
|
||||
assert_(sol.success)
|
||||
|
||||
assert_(sol.x.size < 8)
|
||||
|
||||
sol_test = sol.sol(x_test)
|
||||
assert_allclose(sol_test[0], nonlin_bc_sol(x_test), rtol=1e-5, atol=1e-5)
|
||||
|
||||
f_test = nonlin_bc_fun(x_test, sol_test)
|
||||
r = sol.sol(x_test, 1) - f_test
|
||||
rel_res = r / (1 + np.abs(f_test))
|
||||
norm_res = np.sum(rel_res ** 2, axis=0) ** 0.5
|
||||
|
||||
assert_(np.all(norm_res < 1e-3))
|
||||
assert_allclose(sol.sol(sol.x), sol.y, rtol=1e-10, atol=1e-10)
|
||||
assert_allclose(sol.sol(sol.x, 1), sol.yp, rtol=1e-10, atol=1e-10)
|
||||
|
||||
|
||||
@pytest.mark.thread_unsafe
|
||||
def test_verbose():
|
||||
# Smoke test that checks the printing does something and does not crash
|
||||
x = np.linspace(0, 1, 5)
|
||||
y = np.zeros((2, x.shape[0]))
|
||||
for verbose in [0, 1, 2]:
|
||||
old_stdout = sys.stdout
|
||||
sys.stdout = StringIO()
|
||||
try:
|
||||
sol = solve_bvp(exp_fun, exp_bc, x, y, verbose=verbose)
|
||||
text = sys.stdout.getvalue()
|
||||
finally:
|
||||
sys.stdout = old_stdout
|
||||
|
||||
assert_(sol.success)
|
||||
if verbose == 0:
|
||||
assert_(not text, text)
|
||||
if verbose >= 1:
|
||||
assert_("Solved in" in text, text)
|
||||
if verbose >= 2:
|
||||
assert_("Max residual" in text, text)
|
||||
File diff suppressed because it is too large
Load diff
|
|
@ -1,840 +0,0 @@
|
|||
# Authors: Nils Wagner, Ed Schofield, Pauli Virtanen, John Travers
|
||||
"""
|
||||
Tests for numerical integration.
|
||||
"""
|
||||
import numpy as np
|
||||
from numpy import (arange, zeros, array, dot, sqrt, cos, sin, eye, pi, exp,
|
||||
allclose)
|
||||
|
||||
from numpy.testing import (
|
||||
assert_, assert_array_almost_equal,
|
||||
assert_allclose, assert_array_equal, assert_equal, assert_warns)
|
||||
import pytest
|
||||
from pytest import raises as assert_raises
|
||||
from scipy.integrate import odeint, ode, complex_ode
|
||||
|
||||
#------------------------------------------------------------------------------
|
||||
# Test ODE integrators
|
||||
#------------------------------------------------------------------------------
|
||||
|
||||
|
||||
class TestOdeint:
|
||||
# Check integrate.odeint
|
||||
|
||||
def _do_problem(self, problem):
|
||||
t = arange(0.0, problem.stop_t, 0.05)
|
||||
|
||||
# Basic case
|
||||
z, infodict = odeint(problem.f, problem.z0, t, full_output=True)
|
||||
assert_(problem.verify(z, t))
|
||||
|
||||
# Use tfirst=True
|
||||
z, infodict = odeint(lambda t, y: problem.f(y, t), problem.z0, t,
|
||||
full_output=True, tfirst=True)
|
||||
assert_(problem.verify(z, t))
|
||||
|
||||
if hasattr(problem, 'jac'):
|
||||
# Use Dfun
|
||||
z, infodict = odeint(problem.f, problem.z0, t, Dfun=problem.jac,
|
||||
full_output=True)
|
||||
assert_(problem.verify(z, t))
|
||||
|
||||
# Use Dfun and tfirst=True
|
||||
z, infodict = odeint(lambda t, y: problem.f(y, t), problem.z0, t,
|
||||
Dfun=lambda t, y: problem.jac(y, t),
|
||||
full_output=True, tfirst=True)
|
||||
assert_(problem.verify(z, t))
|
||||
|
||||
def test_odeint(self):
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
if problem.cmplx:
|
||||
continue
|
||||
self._do_problem(problem)
|
||||
|
||||
|
||||
class TestODEClass:
|
||||
|
||||
ode_class = None # Set in subclass.
|
||||
|
||||
def _do_problem(self, problem, integrator, method='adams'):
|
||||
|
||||
# ode has callback arguments in different order than odeint
|
||||
def f(t, z):
|
||||
return problem.f(z, t)
|
||||
jac = None
|
||||
if hasattr(problem, 'jac'):
|
||||
def jac(t, z):
|
||||
return problem.jac(z, t)
|
||||
|
||||
integrator_params = {}
|
||||
if problem.lband is not None or problem.uband is not None:
|
||||
integrator_params['uband'] = problem.uband
|
||||
integrator_params['lband'] = problem.lband
|
||||
|
||||
ig = self.ode_class(f, jac)
|
||||
ig.set_integrator(integrator,
|
||||
atol=problem.atol/10,
|
||||
rtol=problem.rtol/10,
|
||||
method=method,
|
||||
**integrator_params)
|
||||
|
||||
ig.set_initial_value(problem.z0, t=0.0)
|
||||
z = ig.integrate(problem.stop_t)
|
||||
|
||||
assert_array_equal(z, ig.y)
|
||||
assert_(ig.successful(), (problem, method))
|
||||
assert_(ig.get_return_code() > 0, (problem, method))
|
||||
assert_(problem.verify(array([z]), problem.stop_t), (problem, method))
|
||||
|
||||
|
||||
class TestOde(TestODEClass):
|
||||
|
||||
ode_class = ode
|
||||
|
||||
def test_vode(self):
|
||||
# Check the vode solver
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
if problem.cmplx:
|
||||
continue
|
||||
if not problem.stiff:
|
||||
self._do_problem(problem, 'vode', 'adams')
|
||||
self._do_problem(problem, 'vode', 'bdf')
|
||||
|
||||
def test_zvode(self):
|
||||
# Check the zvode solver
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
if not problem.stiff:
|
||||
self._do_problem(problem, 'zvode', 'adams')
|
||||
self._do_problem(problem, 'zvode', 'bdf')
|
||||
|
||||
def test_lsoda(self):
|
||||
# Check the lsoda solver
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
if problem.cmplx:
|
||||
continue
|
||||
self._do_problem(problem, 'lsoda')
|
||||
|
||||
def test_dopri5(self):
|
||||
# Check the dopri5 solver
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
if problem.cmplx:
|
||||
continue
|
||||
if problem.stiff:
|
||||
continue
|
||||
if hasattr(problem, 'jac'):
|
||||
continue
|
||||
self._do_problem(problem, 'dopri5')
|
||||
|
||||
def test_dop853(self):
|
||||
# Check the dop853 solver
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
if problem.cmplx:
|
||||
continue
|
||||
if problem.stiff:
|
||||
continue
|
||||
if hasattr(problem, 'jac'):
|
||||
continue
|
||||
self._do_problem(problem, 'dop853')
|
||||
|
||||
@pytest.mark.thread_unsafe
|
||||
def test_concurrent_fail(self):
|
||||
for sol in ('vode', 'zvode', 'lsoda'):
|
||||
def f(t, y):
|
||||
return 1.0
|
||||
|
||||
r = ode(f).set_integrator(sol)
|
||||
r.set_initial_value(0, 0)
|
||||
|
||||
r2 = ode(f).set_integrator(sol)
|
||||
r2.set_initial_value(0, 0)
|
||||
|
||||
r.integrate(r.t + 0.1)
|
||||
r2.integrate(r2.t + 0.1)
|
||||
|
||||
assert_raises(RuntimeError, r.integrate, r.t + 0.1)
|
||||
|
||||
def test_concurrent_ok(self, num_parallel_threads):
|
||||
def f(t, y):
|
||||
return 1.0
|
||||
|
||||
for k in range(3):
|
||||
for sol in ('vode', 'zvode', 'lsoda', 'dopri5', 'dop853'):
|
||||
if sol in {'vode', 'zvode', 'lsoda'} and num_parallel_threads > 1:
|
||||
continue
|
||||
r = ode(f).set_integrator(sol)
|
||||
r.set_initial_value(0, 0)
|
||||
|
||||
r2 = ode(f).set_integrator(sol)
|
||||
r2.set_initial_value(0, 0)
|
||||
|
||||
r.integrate(r.t + 0.1)
|
||||
r2.integrate(r2.t + 0.1)
|
||||
r2.integrate(r2.t + 0.1)
|
||||
|
||||
assert_allclose(r.y, 0.1)
|
||||
assert_allclose(r2.y, 0.2)
|
||||
|
||||
for sol in ('dopri5', 'dop853'):
|
||||
r = ode(f).set_integrator(sol)
|
||||
r.set_initial_value(0, 0)
|
||||
|
||||
r2 = ode(f).set_integrator(sol)
|
||||
r2.set_initial_value(0, 0)
|
||||
|
||||
r.integrate(r.t + 0.1)
|
||||
r.integrate(r.t + 0.1)
|
||||
r2.integrate(r2.t + 0.1)
|
||||
r.integrate(r.t + 0.1)
|
||||
r2.integrate(r2.t + 0.1)
|
||||
|
||||
assert_allclose(r.y, 0.3)
|
||||
assert_allclose(r2.y, 0.2)
|
||||
|
||||
|
||||
class TestComplexOde(TestODEClass):
|
||||
|
||||
ode_class = complex_ode
|
||||
|
||||
def test_vode(self):
|
||||
# Check the vode solver
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
if not problem.stiff:
|
||||
self._do_problem(problem, 'vode', 'adams')
|
||||
else:
|
||||
self._do_problem(problem, 'vode', 'bdf')
|
||||
|
||||
def test_lsoda(self):
|
||||
|
||||
# Check the lsoda solver
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
self._do_problem(problem, 'lsoda')
|
||||
|
||||
def test_dopri5(self):
|
||||
# Check the dopri5 solver
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
if problem.stiff:
|
||||
continue
|
||||
if hasattr(problem, 'jac'):
|
||||
continue
|
||||
self._do_problem(problem, 'dopri5')
|
||||
|
||||
def test_dop853(self):
|
||||
# Check the dop853 solver
|
||||
for problem_cls in PROBLEMS:
|
||||
problem = problem_cls()
|
||||
if problem.stiff:
|
||||
continue
|
||||
if hasattr(problem, 'jac'):
|
||||
continue
|
||||
self._do_problem(problem, 'dop853')
|
||||
|
||||
|
||||
class TestSolout:
|
||||
# Check integrate.ode correctly handles solout for dopri5 and dop853
|
||||
def _run_solout_test(self, integrator):
|
||||
# Check correct usage of solout
|
||||
ts = []
|
||||
ys = []
|
||||
t0 = 0.0
|
||||
tend = 10.0
|
||||
y0 = [1.0, 2.0]
|
||||
|
||||
def solout(t, y):
|
||||
ts.append(t)
|
||||
ys.append(y.copy())
|
||||
|
||||
def rhs(t, y):
|
||||
return [y[0] + y[1], -y[1]**2]
|
||||
|
||||
ig = ode(rhs).set_integrator(integrator)
|
||||
ig.set_solout(solout)
|
||||
ig.set_initial_value(y0, t0)
|
||||
ret = ig.integrate(tend)
|
||||
assert_array_equal(ys[0], y0)
|
||||
assert_array_equal(ys[-1], ret)
|
||||
assert_equal(ts[0], t0)
|
||||
assert_equal(ts[-1], tend)
|
||||
|
||||
def test_solout(self):
|
||||
for integrator in ('dopri5', 'dop853'):
|
||||
self._run_solout_test(integrator)
|
||||
|
||||
def _run_solout_after_initial_test(self, integrator):
|
||||
# Check if solout works even if it is set after the initial value.
|
||||
ts = []
|
||||
ys = []
|
||||
t0 = 0.0
|
||||
tend = 10.0
|
||||
y0 = [1.0, 2.0]
|
||||
|
||||
def solout(t, y):
|
||||
ts.append(t)
|
||||
ys.append(y.copy())
|
||||
|
||||
def rhs(t, y):
|
||||
return [y[0] + y[1], -y[1]**2]
|
||||
|
||||
ig = ode(rhs).set_integrator(integrator)
|
||||
ig.set_initial_value(y0, t0)
|
||||
ig.set_solout(solout)
|
||||
ret = ig.integrate(tend)
|
||||
assert_array_equal(ys[0], y0)
|
||||
assert_array_equal(ys[-1], ret)
|
||||
assert_equal(ts[0], t0)
|
||||
assert_equal(ts[-1], tend)
|
||||
|
||||
def test_solout_after_initial(self):
|
||||
for integrator in ('dopri5', 'dop853'):
|
||||
self._run_solout_after_initial_test(integrator)
|
||||
|
||||
def _run_solout_break_test(self, integrator):
|
||||
# Check correct usage of stopping via solout
|
||||
ts = []
|
||||
ys = []
|
||||
t0 = 0.0
|
||||
tend = 10.0
|
||||
y0 = [1.0, 2.0]
|
||||
|
||||
def solout(t, y):
|
||||
ts.append(t)
|
||||
ys.append(y.copy())
|
||||
if t > tend/2.0:
|
||||
return -1
|
||||
|
||||
def rhs(t, y):
|
||||
return [y[0] + y[1], -y[1]**2]
|
||||
|
||||
ig = ode(rhs).set_integrator(integrator)
|
||||
ig.set_solout(solout)
|
||||
ig.set_initial_value(y0, t0)
|
||||
ret = ig.integrate(tend)
|
||||
assert_array_equal(ys[0], y0)
|
||||
assert_array_equal(ys[-1], ret)
|
||||
assert_equal(ts[0], t0)
|
||||
assert_(ts[-1] > tend/2.0)
|
||||
assert_(ts[-1] < tend)
|
||||
|
||||
def test_solout_break(self):
|
||||
for integrator in ('dopri5', 'dop853'):
|
||||
self._run_solout_break_test(integrator)
|
||||
|
||||
|
||||
class TestComplexSolout:
|
||||
# Check integrate.ode correctly handles solout for dopri5 and dop853
|
||||
def _run_solout_test(self, integrator):
|
||||
# Check correct usage of solout
|
||||
ts = []
|
||||
ys = []
|
||||
t0 = 0.0
|
||||
tend = 20.0
|
||||
y0 = [0.0]
|
||||
|
||||
def solout(t, y):
|
||||
ts.append(t)
|
||||
ys.append(y.copy())
|
||||
|
||||
def rhs(t, y):
|
||||
return [1.0/(t - 10.0 - 1j)]
|
||||
|
||||
ig = complex_ode(rhs).set_integrator(integrator)
|
||||
ig.set_solout(solout)
|
||||
ig.set_initial_value(y0, t0)
|
||||
ret = ig.integrate(tend)
|
||||
assert_array_equal(ys[0], y0)
|
||||
assert_array_equal(ys[-1], ret)
|
||||
assert_equal(ts[0], t0)
|
||||
assert_equal(ts[-1], tend)
|
||||
|
||||
def test_solout(self):
|
||||
for integrator in ('dopri5', 'dop853'):
|
||||
self._run_solout_test(integrator)
|
||||
|
||||
def _run_solout_break_test(self, integrator):
|
||||
# Check correct usage of stopping via solout
|
||||
ts = []
|
||||
ys = []
|
||||
t0 = 0.0
|
||||
tend = 20.0
|
||||
y0 = [0.0]
|
||||
|
||||
def solout(t, y):
|
||||
ts.append(t)
|
||||
ys.append(y.copy())
|
||||
if t > tend/2.0:
|
||||
return -1
|
||||
|
||||
def rhs(t, y):
|
||||
return [1.0/(t - 10.0 - 1j)]
|
||||
|
||||
ig = complex_ode(rhs).set_integrator(integrator)
|
||||
ig.set_solout(solout)
|
||||
ig.set_initial_value(y0, t0)
|
||||
ret = ig.integrate(tend)
|
||||
assert_array_equal(ys[0], y0)
|
||||
assert_array_equal(ys[-1], ret)
|
||||
assert_equal(ts[0], t0)
|
||||
assert_(ts[-1] > tend/2.0)
|
||||
assert_(ts[-1] < tend)
|
||||
|
||||
def test_solout_break(self):
|
||||
for integrator in ('dopri5', 'dop853'):
|
||||
self._run_solout_break_test(integrator)
|
||||
|
||||
|
||||
#------------------------------------------------------------------------------
|
||||
# Test problems
|
||||
#------------------------------------------------------------------------------
|
||||
|
||||
|
||||
class ODE:
|
||||
"""
|
||||
ODE problem
|
||||
"""
|
||||
stiff = False
|
||||
cmplx = False
|
||||
stop_t = 1
|
||||
z0 = []
|
||||
|
||||
lband = None
|
||||
uband = None
|
||||
|
||||
atol = 1e-6
|
||||
rtol = 1e-5
|
||||
|
||||
|
||||
class SimpleOscillator(ODE):
|
||||
r"""
|
||||
Free vibration of a simple oscillator::
|
||||
m \ddot{u} + k u = 0, u(0) = u_0 \dot{u}(0) \dot{u}_0
|
||||
Solution::
|
||||
u(t) = u_0*cos(sqrt(k/m)*t)+\dot{u}_0*sin(sqrt(k/m)*t)/sqrt(k/m)
|
||||
"""
|
||||
stop_t = 1 + 0.09
|
||||
z0 = array([1.0, 0.1], float)
|
||||
|
||||
k = 4.0
|
||||
m = 1.0
|
||||
|
||||
def f(self, z, t):
|
||||
tmp = zeros((2, 2), float)
|
||||
tmp[0, 1] = 1.0
|
||||
tmp[1, 0] = -self.k / self.m
|
||||
return dot(tmp, z)
|
||||
|
||||
def verify(self, zs, t):
|
||||
omega = sqrt(self.k / self.m)
|
||||
u = self.z0[0]*cos(omega*t) + self.z0[1]*sin(omega*t)/omega
|
||||
return allclose(u, zs[:, 0], atol=self.atol, rtol=self.rtol)
|
||||
|
||||
|
||||
class ComplexExp(ODE):
|
||||
r"""The equation :lm:`\dot u = i u`"""
|
||||
stop_t = 1.23*pi
|
||||
z0 = exp([1j, 2j, 3j, 4j, 5j])
|
||||
cmplx = True
|
||||
|
||||
def f(self, z, t):
|
||||
return 1j*z
|
||||
|
||||
def jac(self, z, t):
|
||||
return 1j*eye(5)
|
||||
|
||||
def verify(self, zs, t):
|
||||
u = self.z0 * exp(1j*t)
|
||||
return allclose(u, zs, atol=self.atol, rtol=self.rtol)
|
||||
|
||||
|
||||
class Pi(ODE):
|
||||
r"""Integrate 1/(t + 1j) from t=-10 to t=10"""
|
||||
stop_t = 20
|
||||
z0 = [0]
|
||||
cmplx = True
|
||||
|
||||
def f(self, z, t):
|
||||
return array([1./(t - 10 + 1j)])
|
||||
|
||||
def verify(self, zs, t):
|
||||
u = -2j * np.arctan(10)
|
||||
return allclose(u, zs[-1, :], atol=self.atol, rtol=self.rtol)
|
||||
|
||||
|
||||
class CoupledDecay(ODE):
|
||||
r"""
|
||||
3 coupled decays suited for banded treatment
|
||||
(banded mode makes it necessary when N>>3)
|
||||
"""
|
||||
|
||||
stiff = True
|
||||
stop_t = 0.5
|
||||
z0 = [5.0, 7.0, 13.0]
|
||||
lband = 1
|
||||
uband = 0
|
||||
|
||||
lmbd = [0.17, 0.23, 0.29] # fictitious decay constants
|
||||
|
||||
def f(self, z, t):
|
||||
lmbd = self.lmbd
|
||||
return np.array([-lmbd[0]*z[0],
|
||||
-lmbd[1]*z[1] + lmbd[0]*z[0],
|
||||
-lmbd[2]*z[2] + lmbd[1]*z[1]])
|
||||
|
||||
def jac(self, z, t):
|
||||
# The full Jacobian is
|
||||
#
|
||||
# [-lmbd[0] 0 0 ]
|
||||
# [ lmbd[0] -lmbd[1] 0 ]
|
||||
# [ 0 lmbd[1] -lmbd[2]]
|
||||
#
|
||||
# The lower and upper bandwidths are lband=1 and uband=0, resp.
|
||||
# The representation of this array in packed format is
|
||||
#
|
||||
# [-lmbd[0] -lmbd[1] -lmbd[2]]
|
||||
# [ lmbd[0] lmbd[1] 0 ]
|
||||
|
||||
lmbd = self.lmbd
|
||||
j = np.zeros((self.lband + self.uband + 1, 3), order='F')
|
||||
|
||||
def set_j(ri, ci, val):
|
||||
j[self.uband + ri - ci, ci] = val
|
||||
set_j(0, 0, -lmbd[0])
|
||||
set_j(1, 0, lmbd[0])
|
||||
set_j(1, 1, -lmbd[1])
|
||||
set_j(2, 1, lmbd[1])
|
||||
set_j(2, 2, -lmbd[2])
|
||||
return j
|
||||
|
||||
def verify(self, zs, t):
|
||||
# Formulae derived by hand
|
||||
lmbd = np.array(self.lmbd)
|
||||
d10 = lmbd[1] - lmbd[0]
|
||||
d21 = lmbd[2] - lmbd[1]
|
||||
d20 = lmbd[2] - lmbd[0]
|
||||
e0 = np.exp(-lmbd[0] * t)
|
||||
e1 = np.exp(-lmbd[1] * t)
|
||||
e2 = np.exp(-lmbd[2] * t)
|
||||
u = np.vstack((
|
||||
self.z0[0] * e0,
|
||||
self.z0[1] * e1 + self.z0[0] * lmbd[0] / d10 * (e0 - e1),
|
||||
self.z0[2] * e2 + self.z0[1] * lmbd[1] / d21 * (e1 - e2) +
|
||||
lmbd[1] * lmbd[0] * self.z0[0] / d10 *
|
||||
(1 / d20 * (e0 - e2) - 1 / d21 * (e1 - e2)))).transpose()
|
||||
return allclose(u, zs, atol=self.atol, rtol=self.rtol)
|
||||
|
||||
|
||||
PROBLEMS = [SimpleOscillator, ComplexExp, Pi, CoupledDecay]
|
||||
|
||||
#------------------------------------------------------------------------------
|
||||
|
||||
|
||||
def f(t, x):
|
||||
dxdt = [x[1], -x[0]]
|
||||
return dxdt
|
||||
|
||||
|
||||
def jac(t, x):
|
||||
j = array([[0.0, 1.0],
|
||||
[-1.0, 0.0]])
|
||||
return j
|
||||
|
||||
|
||||
def f1(t, x, omega):
|
||||
dxdt = [omega*x[1], -omega*x[0]]
|
||||
return dxdt
|
||||
|
||||
|
||||
def jac1(t, x, omega):
|
||||
j = array([[0.0, omega],
|
||||
[-omega, 0.0]])
|
||||
return j
|
||||
|
||||
|
||||
def f2(t, x, omega1, omega2):
|
||||
dxdt = [omega1*x[1], -omega2*x[0]]
|
||||
return dxdt
|
||||
|
||||
|
||||
def jac2(t, x, omega1, omega2):
|
||||
j = array([[0.0, omega1],
|
||||
[-omega2, 0.0]])
|
||||
return j
|
||||
|
||||
|
||||
def fv(t, x, omega):
|
||||
dxdt = [omega[0]*x[1], -omega[1]*x[0]]
|
||||
return dxdt
|
||||
|
||||
|
||||
def jacv(t, x, omega):
|
||||
j = array([[0.0, omega[0]],
|
||||
[-omega[1], 0.0]])
|
||||
return j
|
||||
|
||||
|
||||
class ODECheckParameterUse:
|
||||
"""Call an ode-class solver with several cases of parameter use."""
|
||||
|
||||
# solver_name must be set before tests can be run with this class.
|
||||
|
||||
# Set these in subclasses.
|
||||
solver_name = ''
|
||||
solver_uses_jac = False
|
||||
|
||||
def _get_solver(self, f, jac):
|
||||
solver = ode(f, jac)
|
||||
if self.solver_uses_jac:
|
||||
solver.set_integrator(self.solver_name, atol=1e-9, rtol=1e-7,
|
||||
with_jacobian=self.solver_uses_jac)
|
||||
else:
|
||||
# XXX Shouldn't set_integrator *always* accept the keyword arg
|
||||
# 'with_jacobian', and perhaps raise an exception if it is set
|
||||
# to True if the solver can't actually use it?
|
||||
solver.set_integrator(self.solver_name, atol=1e-9, rtol=1e-7)
|
||||
return solver
|
||||
|
||||
def _check_solver(self, solver):
|
||||
ic = [1.0, 0.0]
|
||||
solver.set_initial_value(ic, 0.0)
|
||||
solver.integrate(pi)
|
||||
assert_array_almost_equal(solver.y, [-1.0, 0.0])
|
||||
|
||||
def test_no_params(self):
|
||||
solver = self._get_solver(f, jac)
|
||||
self._check_solver(solver)
|
||||
|
||||
def test_one_scalar_param(self):
|
||||
solver = self._get_solver(f1, jac1)
|
||||
omega = 1.0
|
||||
solver.set_f_params(omega)
|
||||
if self.solver_uses_jac:
|
||||
solver.set_jac_params(omega)
|
||||
self._check_solver(solver)
|
||||
|
||||
def test_two_scalar_params(self):
|
||||
solver = self._get_solver(f2, jac2)
|
||||
omega1 = 1.0
|
||||
omega2 = 1.0
|
||||
solver.set_f_params(omega1, omega2)
|
||||
if self.solver_uses_jac:
|
||||
solver.set_jac_params(omega1, omega2)
|
||||
self._check_solver(solver)
|
||||
|
||||
def test_vector_param(self):
|
||||
solver = self._get_solver(fv, jacv)
|
||||
omega = [1.0, 1.0]
|
||||
solver.set_f_params(omega)
|
||||
if self.solver_uses_jac:
|
||||
solver.set_jac_params(omega)
|
||||
self._check_solver(solver)
|
||||
|
||||
@pytest.mark.thread_unsafe
|
||||
def test_warns_on_failure(self):
|
||||
# Set nsteps small to ensure failure
|
||||
solver = self._get_solver(f, jac)
|
||||
solver.set_integrator(self.solver_name, nsteps=1)
|
||||
ic = [1.0, 0.0]
|
||||
solver.set_initial_value(ic, 0.0)
|
||||
assert_warns(UserWarning, solver.integrate, pi)
|
||||
|
||||
|
||||
class TestDOPRI5CheckParameterUse(ODECheckParameterUse):
|
||||
solver_name = 'dopri5'
|
||||
solver_uses_jac = False
|
||||
|
||||
|
||||
class TestDOP853CheckParameterUse(ODECheckParameterUse):
|
||||
solver_name = 'dop853'
|
||||
solver_uses_jac = False
|
||||
|
||||
|
||||
class TestVODECheckParameterUse(ODECheckParameterUse):
|
||||
solver_name = 'vode'
|
||||
solver_uses_jac = True
|
||||
|
||||
|
||||
class TestZVODECheckParameterUse(ODECheckParameterUse):
|
||||
solver_name = 'zvode'
|
||||
solver_uses_jac = True
|
||||
|
||||
|
||||
class TestLSODACheckParameterUse(ODECheckParameterUse):
|
||||
solver_name = 'lsoda'
|
||||
solver_uses_jac = True
|
||||
|
||||
|
||||
def test_odeint_trivial_time():
|
||||
# Test that odeint succeeds when given a single time point
|
||||
# and full_output=True. This is a regression test for gh-4282.
|
||||
y0 = 1
|
||||
t = [0]
|
||||
y, info = odeint(lambda y, t: -y, y0, t, full_output=True)
|
||||
assert_array_equal(y, np.array([[y0]]))
|
||||
|
||||
|
||||
def test_odeint_banded_jacobian():
|
||||
# Test the use of the `Dfun`, `ml` and `mu` options of odeint.
|
||||
|
||||
def func(y, t, c):
|
||||
return c.dot(y)
|
||||
|
||||
def jac(y, t, c):
|
||||
return c
|
||||
|
||||
def jac_transpose(y, t, c):
|
||||
return c.T.copy(order='C')
|
||||
|
||||
def bjac_rows(y, t, c):
|
||||
jac = np.vstack((np.r_[0, np.diag(c, 1)],
|
||||
np.diag(c),
|
||||
np.r_[np.diag(c, -1), 0],
|
||||
np.r_[np.diag(c, -2), 0, 0]))
|
||||
return jac
|
||||
|
||||
def bjac_cols(y, t, c):
|
||||
return bjac_rows(y, t, c).T.copy(order='C')
|
||||
|
||||
c = array([[-205, 0.01, 0.00, 0.0],
|
||||
[0.1, -2.50, 0.02, 0.0],
|
||||
[1e-3, 0.01, -2.0, 0.01],
|
||||
[0.00, 0.00, 0.1, -1.0]])
|
||||
|
||||
y0 = np.ones(4)
|
||||
t = np.array([0, 5, 10, 100])
|
||||
|
||||
# Use the full Jacobian.
|
||||
sol1, info1 = odeint(func, y0, t, args=(c,), full_output=True,
|
||||
atol=1e-13, rtol=1e-11, mxstep=10000,
|
||||
Dfun=jac)
|
||||
|
||||
# Use the transposed full Jacobian, with col_deriv=True.
|
||||
sol2, info2 = odeint(func, y0, t, args=(c,), full_output=True,
|
||||
atol=1e-13, rtol=1e-11, mxstep=10000,
|
||||
Dfun=jac_transpose, col_deriv=True)
|
||||
|
||||
# Use the banded Jacobian.
|
||||
sol3, info3 = odeint(func, y0, t, args=(c,), full_output=True,
|
||||
atol=1e-13, rtol=1e-11, mxstep=10000,
|
||||
Dfun=bjac_rows, ml=2, mu=1)
|
||||
|
||||
# Use the transposed banded Jacobian, with col_deriv=True.
|
||||
sol4, info4 = odeint(func, y0, t, args=(c,), full_output=True,
|
||||
atol=1e-13, rtol=1e-11, mxstep=10000,
|
||||
Dfun=bjac_cols, ml=2, mu=1, col_deriv=True)
|
||||
|
||||
assert_allclose(sol1, sol2, err_msg="sol1 != sol2")
|
||||
assert_allclose(sol1, sol3, atol=1e-12, err_msg="sol1 != sol3")
|
||||
assert_allclose(sol3, sol4, err_msg="sol3 != sol4")
|
||||
|
||||
# Verify that the number of jacobian evaluations was the same for the
|
||||
# calls of odeint with a full jacobian and with a banded jacobian. This is
|
||||
# a regression test--there was a bug in the handling of banded jacobians
|
||||
# that resulted in an incorrect jacobian matrix being passed to the LSODA
|
||||
# code. That would cause errors or excessive jacobian evaluations.
|
||||
assert_array_equal(info1['nje'], info2['nje'])
|
||||
assert_array_equal(info3['nje'], info4['nje'])
|
||||
|
||||
# Test the use of tfirst
|
||||
sol1ty, info1ty = odeint(lambda t, y, c: func(y, t, c), y0, t, args=(c,),
|
||||
full_output=True, atol=1e-13, rtol=1e-11,
|
||||
mxstep=10000,
|
||||
Dfun=lambda t, y, c: jac(y, t, c), tfirst=True)
|
||||
# The code should execute the exact same sequence of floating point
|
||||
# calculations, so these should be exactly equal. We'll be safe and use
|
||||
# a small tolerance.
|
||||
assert_allclose(sol1, sol1ty, rtol=1e-12, err_msg="sol1 != sol1ty")
|
||||
|
||||
|
||||
def test_odeint_errors():
|
||||
def sys1d(x, t):
|
||||
return -100*x
|
||||
|
||||
def bad1(x, t):
|
||||
return 1.0/0
|
||||
|
||||
def bad2(x, t):
|
||||
return "foo"
|
||||
|
||||
def bad_jac1(x, t):
|
||||
return 1.0/0
|
||||
|
||||
def bad_jac2(x, t):
|
||||
return [["foo"]]
|
||||
|
||||
def sys2d(x, t):
|
||||
return [-100*x[0], -0.1*x[1]]
|
||||
|
||||
def sys2d_bad_jac(x, t):
|
||||
return [[1.0/0, 0], [0, -0.1]]
|
||||
|
||||
assert_raises(ZeroDivisionError, odeint, bad1, 1.0, [0, 1])
|
||||
assert_raises(ValueError, odeint, bad2, 1.0, [0, 1])
|
||||
|
||||
assert_raises(ZeroDivisionError, odeint, sys1d, 1.0, [0, 1], Dfun=bad_jac1)
|
||||
assert_raises(ValueError, odeint, sys1d, 1.0, [0, 1], Dfun=bad_jac2)
|
||||
|
||||
assert_raises(ZeroDivisionError, odeint, sys2d, [1.0, 1.0], [0, 1],
|
||||
Dfun=sys2d_bad_jac)
|
||||
|
||||
|
||||
def test_odeint_bad_shapes():
|
||||
# Tests of some errors that can occur with odeint.
|
||||
|
||||
def badrhs(x, t):
|
||||
return [1, -1]
|
||||
|
||||
def sys1(x, t):
|
||||
return -100*x
|
||||
|
||||
def badjac(x, t):
|
||||
return [[0, 0, 0]]
|
||||
|
||||
# y0 must be at most 1-d.
|
||||
bad_y0 = [[0, 0], [0, 0]]
|
||||
assert_raises(ValueError, odeint, sys1, bad_y0, [0, 1])
|
||||
|
||||
# t must be at most 1-d.
|
||||
bad_t = [[0, 1], [2, 3]]
|
||||
assert_raises(ValueError, odeint, sys1, [10.0], bad_t)
|
||||
|
||||
# y0 is 10, but badrhs(x, t) returns [1, -1].
|
||||
assert_raises(RuntimeError, odeint, badrhs, 10, [0, 1])
|
||||
|
||||
# shape of array returned by badjac(x, t) is not correct.
|
||||
assert_raises(RuntimeError, odeint, sys1, [10, 10], [0, 1], Dfun=badjac)
|
||||
|
||||
|
||||
def test_repeated_t_values():
|
||||
"""Regression test for gh-8217."""
|
||||
|
||||
def func(x, t):
|
||||
return -0.25*x
|
||||
|
||||
t = np.zeros(10)
|
||||
sol = odeint(func, [1.], t)
|
||||
assert_array_equal(sol, np.ones((len(t), 1)))
|
||||
|
||||
tau = 4*np.log(2)
|
||||
t = [0]*9 + [tau, 2*tau, 2*tau, 3*tau]
|
||||
sol = odeint(func, [1, 2], t, rtol=1e-12, atol=1e-12)
|
||||
expected_sol = np.array([[1.0, 2.0]]*9 +
|
||||
[[0.5, 1.0],
|
||||
[0.25, 0.5],
|
||||
[0.25, 0.5],
|
||||
[0.125, 0.25]])
|
||||
assert_allclose(sol, expected_sol)
|
||||
|
||||
# Edge case: empty t sequence.
|
||||
sol = odeint(func, [1.], [])
|
||||
assert_array_equal(sol, np.array([], dtype=np.float64).reshape((0, 1)))
|
||||
|
||||
# t values are not monotonic.
|
||||
assert_raises(ValueError, odeint, func, [1.], [0, 1, 0.5, 0])
|
||||
assert_raises(ValueError, odeint, func, [1, 2, 3], [0, -1, -2, 3])
|
||||
|
|
@ -1,74 +0,0 @@
|
|||
import numpy as np
|
||||
from numpy.testing import assert_equal, assert_allclose
|
||||
from scipy.integrate import odeint
|
||||
import scipy.integrate._test_odeint_banded as banded5x5
|
||||
|
||||
|
||||
def rhs(y, t):
|
||||
dydt = np.zeros_like(y)
|
||||
banded5x5.banded5x5(t, y, dydt)
|
||||
return dydt
|
||||
|
||||
|
||||
def jac(y, t):
|
||||
n = len(y)
|
||||
jac = np.zeros((n, n), order='F')
|
||||
banded5x5.banded5x5_jac(t, y, 1, 1, jac)
|
||||
return jac
|
||||
|
||||
|
||||
def bjac(y, t):
|
||||
n = len(y)
|
||||
bjac = np.zeros((4, n), order='F')
|
||||
banded5x5.banded5x5_bjac(t, y, 1, 1, bjac)
|
||||
return bjac
|
||||
|
||||
|
||||
JACTYPE_FULL = 1
|
||||
JACTYPE_BANDED = 4
|
||||
|
||||
|
||||
def check_odeint(jactype):
|
||||
if jactype == JACTYPE_FULL:
|
||||
ml = None
|
||||
mu = None
|
||||
jacobian = jac
|
||||
elif jactype == JACTYPE_BANDED:
|
||||
ml = 2
|
||||
mu = 1
|
||||
jacobian = bjac
|
||||
else:
|
||||
raise ValueError(f"invalid jactype: {jactype!r}")
|
||||
|
||||
y0 = np.arange(1.0, 6.0)
|
||||
# These tolerances must match the tolerances used in banded5x5.f.
|
||||
rtol = 1e-11
|
||||
atol = 1e-13
|
||||
dt = 0.125
|
||||
nsteps = 64
|
||||
t = dt * np.arange(nsteps+1)
|
||||
|
||||
sol, info = odeint(rhs, y0, t,
|
||||
Dfun=jacobian, ml=ml, mu=mu,
|
||||
atol=atol, rtol=rtol, full_output=True)
|
||||
yfinal = sol[-1]
|
||||
odeint_nst = info['nst'][-1]
|
||||
odeint_nfe = info['nfe'][-1]
|
||||
odeint_nje = info['nje'][-1]
|
||||
|
||||
y1 = y0.copy()
|
||||
# Pure Fortran solution. y1 is modified in-place.
|
||||
nst, nfe, nje = banded5x5.banded5x5_solve(y1, nsteps, dt, jactype)
|
||||
|
||||
# It is likely that yfinal and y1 are *exactly* the same, but
|
||||
# we'll be cautious and use assert_allclose.
|
||||
assert_allclose(yfinal, y1, rtol=1e-12)
|
||||
assert_equal((odeint_nst, odeint_nfe, odeint_nje), (nst, nfe, nje))
|
||||
|
||||
|
||||
def test_odeint_full_jac():
|
||||
check_odeint(JACTYPE_FULL)
|
||||
|
||||
|
||||
def test_odeint_banded_jac():
|
||||
check_odeint(JACTYPE_BANDED)
|
||||
|
|
@ -1,680 +0,0 @@
|
|||
import sys
|
||||
import math
|
||||
import numpy as np
|
||||
from numpy import sqrt, cos, sin, arctan, exp, log, pi
|
||||
from numpy.testing import (assert_,
|
||||
assert_allclose, assert_array_less, assert_almost_equal)
|
||||
import pytest
|
||||
|
||||
from scipy.integrate import quad, dblquad, tplquad, nquad
|
||||
from scipy.special import erf, erfc
|
||||
from scipy._lib._ccallback import LowLevelCallable
|
||||
|
||||
import ctypes
|
||||
import ctypes.util
|
||||
from scipy._lib._ccallback_c import sine_ctypes
|
||||
|
||||
import scipy.integrate._test_multivariate as clib_test
|
||||
|
||||
|
||||
def assert_quad(value_and_err, tabled_value, error_tolerance=1.5e-8):
|
||||
value, err = value_and_err
|
||||
assert_allclose(value, tabled_value, atol=err, rtol=0)
|
||||
if error_tolerance is not None:
|
||||
assert_array_less(err, error_tolerance)
|
||||
|
||||
|
||||
def get_clib_test_routine(name, restype, *argtypes):
|
||||
ptr = getattr(clib_test, name)
|
||||
return ctypes.cast(ptr, ctypes.CFUNCTYPE(restype, *argtypes))
|
||||
|
||||
|
||||
class TestCtypesQuad:
|
||||
def setup_method(self):
|
||||
if sys.platform == 'win32':
|
||||
files = ['api-ms-win-crt-math-l1-1-0.dll']
|
||||
elif sys.platform == 'darwin':
|
||||
files = ['libm.dylib']
|
||||
else:
|
||||
files = ['libm.so', 'libm.so.6']
|
||||
|
||||
for file in files:
|
||||
try:
|
||||
self.lib = ctypes.CDLL(file)
|
||||
break
|
||||
except OSError:
|
||||
pass
|
||||
else:
|
||||
# This test doesn't work on some Linux platforms (Fedora for
|
||||
# example) that put an ld script in libm.so - see gh-5370
|
||||
pytest.skip("Ctypes can't import libm.so")
|
||||
|
||||
restype = ctypes.c_double
|
||||
argtypes = (ctypes.c_double,)
|
||||
for name in ['sin', 'cos', 'tan']:
|
||||
func = getattr(self.lib, name)
|
||||
func.restype = restype
|
||||
func.argtypes = argtypes
|
||||
|
||||
def test_typical(self):
|
||||
assert_quad(quad(self.lib.sin, 0, 5), quad(math.sin, 0, 5)[0])
|
||||
assert_quad(quad(self.lib.cos, 0, 5), quad(math.cos, 0, 5)[0])
|
||||
assert_quad(quad(self.lib.tan, 0, 1), quad(math.tan, 0, 1)[0])
|
||||
|
||||
def test_ctypes_sine(self):
|
||||
quad(LowLevelCallable(sine_ctypes), 0, 1)
|
||||
|
||||
def test_ctypes_variants(self):
|
||||
sin_0 = get_clib_test_routine('_sin_0', ctypes.c_double,
|
||||
ctypes.c_double, ctypes.c_void_p)
|
||||
|
||||
sin_1 = get_clib_test_routine('_sin_1', ctypes.c_double,
|
||||
ctypes.c_int, ctypes.POINTER(ctypes.c_double),
|
||||
ctypes.c_void_p)
|
||||
|
||||
sin_2 = get_clib_test_routine('_sin_2', ctypes.c_double,
|
||||
ctypes.c_double)
|
||||
|
||||
sin_3 = get_clib_test_routine('_sin_3', ctypes.c_double,
|
||||
ctypes.c_int, ctypes.POINTER(ctypes.c_double))
|
||||
|
||||
sin_4 = get_clib_test_routine('_sin_3', ctypes.c_double,
|
||||
ctypes.c_int, ctypes.c_double)
|
||||
|
||||
all_sigs = [sin_0, sin_1, sin_2, sin_3, sin_4]
|
||||
legacy_sigs = [sin_2, sin_4]
|
||||
legacy_only_sigs = [sin_4]
|
||||
|
||||
# LowLevelCallables work for new signatures
|
||||
for j, func in enumerate(all_sigs):
|
||||
callback = LowLevelCallable(func)
|
||||
if func in legacy_only_sigs:
|
||||
pytest.raises(ValueError, quad, callback, 0, pi)
|
||||
else:
|
||||
assert_allclose(quad(callback, 0, pi)[0], 2.0)
|
||||
|
||||
# Plain ctypes items work only for legacy signatures
|
||||
for j, func in enumerate(legacy_sigs):
|
||||
if func in legacy_sigs:
|
||||
assert_allclose(quad(func, 0, pi)[0], 2.0)
|
||||
else:
|
||||
pytest.raises(ValueError, quad, func, 0, pi)
|
||||
|
||||
|
||||
class TestMultivariateCtypesQuad:
|
||||
def setup_method(self):
|
||||
restype = ctypes.c_double
|
||||
argtypes = (ctypes.c_int, ctypes.c_double)
|
||||
for name in ['_multivariate_typical', '_multivariate_indefinite',
|
||||
'_multivariate_sin']:
|
||||
func = get_clib_test_routine(name, restype, *argtypes)
|
||||
setattr(self, name, func)
|
||||
|
||||
def test_typical(self):
|
||||
# 1) Typical function with two extra arguments:
|
||||
assert_quad(quad(self._multivariate_typical, 0, pi, (2, 1.8)),
|
||||
0.30614353532540296487)
|
||||
|
||||
def test_indefinite(self):
|
||||
# 2) Infinite integration limits --- Euler's constant
|
||||
assert_quad(quad(self._multivariate_indefinite, 0, np.inf),
|
||||
0.577215664901532860606512)
|
||||
|
||||
def test_threadsafety(self):
|
||||
# Ensure multivariate ctypes are threadsafe
|
||||
def threadsafety(y):
|
||||
return y + quad(self._multivariate_sin, 0, 1)[0]
|
||||
assert_quad(quad(threadsafety, 0, 1), 0.9596976941318602)
|
||||
|
||||
|
||||
class TestQuad:
|
||||
def test_typical(self):
|
||||
# 1) Typical function with two extra arguments:
|
||||
def myfunc(x, n, z): # Bessel function integrand
|
||||
return cos(n*x-z*sin(x))/pi
|
||||
assert_quad(quad(myfunc, 0, pi, (2, 1.8)), 0.30614353532540296487)
|
||||
|
||||
def test_indefinite(self):
|
||||
# 2) Infinite integration limits --- Euler's constant
|
||||
def myfunc(x): # Euler's constant integrand
|
||||
return -exp(-x)*log(x)
|
||||
assert_quad(quad(myfunc, 0, np.inf), 0.577215664901532860606512)
|
||||
|
||||
def test_singular(self):
|
||||
# 3) Singular points in region of integration.
|
||||
def myfunc(x):
|
||||
if 0 < x < 2.5:
|
||||
return sin(x)
|
||||
elif 2.5 <= x <= 5.0:
|
||||
return exp(-x)
|
||||
else:
|
||||
return 0.0
|
||||
|
||||
assert_quad(quad(myfunc, 0, 10, points=[2.5, 5.0]),
|
||||
1 - cos(2.5) + exp(-2.5) - exp(-5.0))
|
||||
|
||||
def test_sine_weighted_finite(self):
|
||||
# 4) Sine weighted integral (finite limits)
|
||||
def myfunc(x, a):
|
||||
return exp(a*(x-1))
|
||||
|
||||
ome = 2.0**3.4
|
||||
assert_quad(quad(myfunc, 0, 1, args=20, weight='sin', wvar=ome),
|
||||
(20*sin(ome)-ome*cos(ome)+ome*exp(-20))/(20**2 + ome**2))
|
||||
|
||||
def test_sine_weighted_infinite(self):
|
||||
# 5) Sine weighted integral (infinite limits)
|
||||
def myfunc(x, a):
|
||||
return exp(-x*a)
|
||||
|
||||
a = 4.0
|
||||
ome = 3.0
|
||||
assert_quad(quad(myfunc, 0, np.inf, args=a, weight='sin', wvar=ome),
|
||||
ome/(a**2 + ome**2))
|
||||
|
||||
def test_cosine_weighted_infinite(self):
|
||||
# 6) Cosine weighted integral (negative infinite limits)
|
||||
def myfunc(x, a):
|
||||
return exp(x*a)
|
||||
|
||||
a = 2.5
|
||||
ome = 2.3
|
||||
assert_quad(quad(myfunc, -np.inf, 0, args=a, weight='cos', wvar=ome),
|
||||
a/(a**2 + ome**2))
|
||||
|
||||
def test_algebraic_log_weight(self):
|
||||
# 6) Algebraic-logarithmic weight.
|
||||
def myfunc(x, a):
|
||||
return 1/(1+x+2**(-a))
|
||||
|
||||
a = 1.5
|
||||
assert_quad(quad(myfunc, -1, 1, args=a, weight='alg',
|
||||
wvar=(-0.5, -0.5)),
|
||||
pi/sqrt((1+2**(-a))**2 - 1))
|
||||
|
||||
def test_cauchypv_weight(self):
|
||||
# 7) Cauchy prinicpal value weighting w(x) = 1/(x-c)
|
||||
def myfunc(x, a):
|
||||
return 2.0**(-a)/((x-1)**2+4.0**(-a))
|
||||
|
||||
a = 0.4
|
||||
tabledValue = ((2.0**(-0.4)*log(1.5) -
|
||||
2.0**(-1.4)*log((4.0**(-a)+16) / (4.0**(-a)+1)) -
|
||||
arctan(2.0**(a+2)) -
|
||||
arctan(2.0**a)) /
|
||||
(4.0**(-a) + 1))
|
||||
assert_quad(quad(myfunc, 0, 5, args=0.4, weight='cauchy', wvar=2.0),
|
||||
tabledValue, error_tolerance=1.9e-8)
|
||||
|
||||
def test_b_less_than_a(self):
|
||||
def f(x, p, q):
|
||||
return p * np.exp(-q*x)
|
||||
|
||||
val_1, err_1 = quad(f, 0, np.inf, args=(2, 3))
|
||||
val_2, err_2 = quad(f, np.inf, 0, args=(2, 3))
|
||||
assert_allclose(val_1, -val_2, atol=max(err_1, err_2))
|
||||
|
||||
def test_b_less_than_a_2(self):
|
||||
def f(x, s):
|
||||
return np.exp(-x**2 / 2 / s) / np.sqrt(2.*s)
|
||||
|
||||
val_1, err_1 = quad(f, -np.inf, np.inf, args=(2,))
|
||||
val_2, err_2 = quad(f, np.inf, -np.inf, args=(2,))
|
||||
assert_allclose(val_1, -val_2, atol=max(err_1, err_2))
|
||||
|
||||
def test_b_less_than_a_3(self):
|
||||
def f(x):
|
||||
return 1.0
|
||||
|
||||
val_1, err_1 = quad(f, 0, 1, weight='alg', wvar=(0, 0))
|
||||
val_2, err_2 = quad(f, 1, 0, weight='alg', wvar=(0, 0))
|
||||
assert_allclose(val_1, -val_2, atol=max(err_1, err_2))
|
||||
|
||||
def test_b_less_than_a_full_output(self):
|
||||
def f(x):
|
||||
return 1.0
|
||||
|
||||
res_1 = quad(f, 0, 1, weight='alg', wvar=(0, 0), full_output=True)
|
||||
res_2 = quad(f, 1, 0, weight='alg', wvar=(0, 0), full_output=True)
|
||||
err = max(res_1[1], res_2[1])
|
||||
assert_allclose(res_1[0], -res_2[0], atol=err)
|
||||
|
||||
def test_double_integral(self):
|
||||
# 8) Double Integral test
|
||||
def simpfunc(y, x): # Note order of arguments.
|
||||
return x+y
|
||||
|
||||
a, b = 1.0, 2.0
|
||||
assert_quad(dblquad(simpfunc, a, b, lambda x: x, lambda x: 2*x),
|
||||
5/6.0 * (b**3.0-a**3.0))
|
||||
|
||||
def test_double_integral2(self):
|
||||
def func(x0, x1, t0, t1):
|
||||
return x0 + x1 + t0 + t1
|
||||
def g(x):
|
||||
return x
|
||||
def h(x):
|
||||
return 2 * x
|
||||
args = 1, 2
|
||||
assert_quad(dblquad(func, 1, 2, g, h, args=args),35./6 + 9*.5)
|
||||
|
||||
def test_double_integral3(self):
|
||||
def func(x0, x1):
|
||||
return x0 + x1 + 1 + 2
|
||||
assert_quad(dblquad(func, 1, 2, 1, 2),6.)
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"x_lower, x_upper, y_lower, y_upper, expected",
|
||||
[
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [-inf, 0] for all n.
|
||||
(-np.inf, 0, -np.inf, 0, np.pi / 4),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [-inf, -1] for each n (one at a time).
|
||||
(-np.inf, -1, -np.inf, 0, np.pi / 4 * erfc(1)),
|
||||
(-np.inf, 0, -np.inf, -1, np.pi / 4 * erfc(1)),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [-inf, -1] for all n.
|
||||
(-np.inf, -1, -np.inf, -1, np.pi / 4 * (erfc(1) ** 2)),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [-inf, 1] for each n (one at a time).
|
||||
(-np.inf, 1, -np.inf, 0, np.pi / 4 * (erf(1) + 1)),
|
||||
(-np.inf, 0, -np.inf, 1, np.pi / 4 * (erf(1) + 1)),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [-inf, 1] for all n.
|
||||
(-np.inf, 1, -np.inf, 1, np.pi / 4 * ((erf(1) + 1) ** 2)),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain Dx = [-inf, -1] and Dy = [-inf, 1].
|
||||
(-np.inf, -1, -np.inf, 1, np.pi / 4 * ((erf(1) + 1) * erfc(1))),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain Dx = [-inf, 1] and Dy = [-inf, -1].
|
||||
(-np.inf, 1, -np.inf, -1, np.pi / 4 * ((erf(1) + 1) * erfc(1))),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [0, inf] for all n.
|
||||
(0, np.inf, 0, np.inf, np.pi / 4),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [1, inf] for each n (one at a time).
|
||||
(1, np.inf, 0, np.inf, np.pi / 4 * erfc(1)),
|
||||
(0, np.inf, 1, np.inf, np.pi / 4 * erfc(1)),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [1, inf] for all n.
|
||||
(1, np.inf, 1, np.inf, np.pi / 4 * (erfc(1) ** 2)),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [-1, inf] for each n (one at a time).
|
||||
(-1, np.inf, 0, np.inf, np.pi / 4 * (erf(1) + 1)),
|
||||
(0, np.inf, -1, np.inf, np.pi / 4 * (erf(1) + 1)),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [-1, inf] for all n.
|
||||
(-1, np.inf, -1, np.inf, np.pi / 4 * ((erf(1) + 1) ** 2)),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain Dx = [-1, inf] and Dy = [1, inf].
|
||||
(-1, np.inf, 1, np.inf, np.pi / 4 * ((erf(1) + 1) * erfc(1))),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain Dx = [1, inf] and Dy = [-1, inf].
|
||||
(1, np.inf, -1, np.inf, np.pi / 4 * ((erf(1) + 1) * erfc(1))),
|
||||
# Multiple integration of a function in n = 2 variables: f(x, y, z)
|
||||
# over domain D = [-inf, inf] for all n.
|
||||
(-np.inf, np.inf, -np.inf, np.inf, np.pi)
|
||||
]
|
||||
)
|
||||
def test_double_integral_improper(
|
||||
self, x_lower, x_upper, y_lower, y_upper, expected
|
||||
):
|
||||
# The Gaussian Integral.
|
||||
def f(x, y):
|
||||
return np.exp(-x ** 2 - y ** 2)
|
||||
|
||||
assert_quad(
|
||||
dblquad(f, x_lower, x_upper, y_lower, y_upper),
|
||||
expected,
|
||||
error_tolerance=3e-8
|
||||
)
|
||||
|
||||
def test_triple_integral(self):
|
||||
# 9) Triple Integral test
|
||||
def simpfunc(z, y, x, t): # Note order of arguments.
|
||||
return (x+y+z)*t
|
||||
|
||||
a, b = 1.0, 2.0
|
||||
assert_quad(tplquad(simpfunc, a, b,
|
||||
lambda x: x, lambda x: 2*x,
|
||||
lambda x, y: x - y, lambda x, y: x + y,
|
||||
(2.,)),
|
||||
2*8/3.0 * (b**4.0 - a**4.0))
|
||||
|
||||
@pytest.mark.xslow
|
||||
@pytest.mark.parametrize(
|
||||
"x_lower, x_upper, y_lower, y_upper, z_lower, z_upper, expected",
|
||||
[
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-inf, 0] for all n.
|
||||
(-np.inf, 0, -np.inf, 0, -np.inf, 0, (np.pi ** (3 / 2)) / 8),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-inf, -1] for each n (one at a time).
|
||||
(-np.inf, -1, -np.inf, 0, -np.inf, 0,
|
||||
(np.pi ** (3 / 2)) / 8 * erfc(1)),
|
||||
(-np.inf, 0, -np.inf, -1, -np.inf, 0,
|
||||
(np.pi ** (3 / 2)) / 8 * erfc(1)),
|
||||
(-np.inf, 0, -np.inf, 0, -np.inf, -1,
|
||||
(np.pi ** (3 / 2)) / 8 * erfc(1)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-inf, -1] for each n (two at a time).
|
||||
(-np.inf, -1, -np.inf, -1, -np.inf, 0,
|
||||
(np.pi ** (3 / 2)) / 8 * (erfc(1) ** 2)),
|
||||
(-np.inf, -1, -np.inf, 0, -np.inf, -1,
|
||||
(np.pi ** (3 / 2)) / 8 * (erfc(1) ** 2)),
|
||||
(-np.inf, 0, -np.inf, -1, -np.inf, -1,
|
||||
(np.pi ** (3 / 2)) / 8 * (erfc(1) ** 2)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-inf, -1] for all n.
|
||||
(-np.inf, -1, -np.inf, -1, -np.inf, -1,
|
||||
(np.pi ** (3 / 2)) / 8 * (erfc(1) ** 3)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = [-inf, -1] and Dy = Dz = [-inf, 1].
|
||||
(-np.inf, -1, -np.inf, 1, -np.inf, 1,
|
||||
(np.pi ** (3 / 2)) / 8 * (((erf(1) + 1) ** 2) * erfc(1))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = Dy = [-inf, -1] and Dz = [-inf, 1].
|
||||
(-np.inf, -1, -np.inf, -1, -np.inf, 1,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) * (erfc(1) ** 2))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = Dz = [-inf, -1] and Dy = [-inf, 1].
|
||||
(-np.inf, -1, -np.inf, 1, -np.inf, -1,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) * (erfc(1) ** 2))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = [-inf, 1] and Dy = Dz = [-inf, -1].
|
||||
(-np.inf, 1, -np.inf, -1, -np.inf, -1,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) * (erfc(1) ** 2))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = Dy = [-inf, 1] and Dz = [-inf, -1].
|
||||
(-np.inf, 1, -np.inf, 1, -np.inf, -1,
|
||||
(np.pi ** (3 / 2)) / 8 * (((erf(1) + 1) ** 2) * erfc(1))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = Dz = [-inf, 1] and Dy = [-inf, -1].
|
||||
(-np.inf, 1, -np.inf, -1, -np.inf, 1,
|
||||
(np.pi ** (3 / 2)) / 8 * (((erf(1) + 1) ** 2) * erfc(1))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-inf, 1] for each n (one at a time).
|
||||
(-np.inf, 1, -np.inf, 0, -np.inf, 0,
|
||||
(np.pi ** (3 / 2)) / 8 * (erf(1) + 1)),
|
||||
(-np.inf, 0, -np.inf, 1, -np.inf, 0,
|
||||
(np.pi ** (3 / 2)) / 8 * (erf(1) + 1)),
|
||||
(-np.inf, 0, -np.inf, 0, -np.inf, 1,
|
||||
(np.pi ** (3 / 2)) / 8 * (erf(1) + 1)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-inf, 1] for each n (two at a time).
|
||||
(-np.inf, 1, -np.inf, 1, -np.inf, 0,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) ** 2)),
|
||||
(-np.inf, 1, -np.inf, 0, -np.inf, 1,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) ** 2)),
|
||||
(-np.inf, 0, -np.inf, 1, -np.inf, 1,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) ** 2)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-inf, 1] for all n.
|
||||
(-np.inf, 1, -np.inf, 1, -np.inf, 1,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) ** 3)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [0, inf] for all n.
|
||||
(0, np.inf, 0, np.inf, 0, np.inf, (np.pi ** (3 / 2)) / 8),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [1, inf] for each n (one at a time).
|
||||
(1, np.inf, 0, np.inf, 0, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * erfc(1)),
|
||||
(0, np.inf, 1, np.inf, 0, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * erfc(1)),
|
||||
(0, np.inf, 0, np.inf, 1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * erfc(1)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [1, inf] for each n (two at a time).
|
||||
(1, np.inf, 1, np.inf, 0, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (erfc(1) ** 2)),
|
||||
(1, np.inf, 0, np.inf, 1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (erfc(1) ** 2)),
|
||||
(0, np.inf, 1, np.inf, 1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (erfc(1) ** 2)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [1, inf] for all n.
|
||||
(1, np.inf, 1, np.inf, 1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (erfc(1) ** 3)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-1, inf] for each n (one at a time).
|
||||
(-1, np.inf, 0, np.inf, 0, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (erf(1) + 1)),
|
||||
(0, np.inf, -1, np.inf, 0, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (erf(1) + 1)),
|
||||
(0, np.inf, 0, np.inf, -1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (erf(1) + 1)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-1, inf] for each n (two at a time).
|
||||
(-1, np.inf, -1, np.inf, 0, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) ** 2)),
|
||||
(-1, np.inf, 0, np.inf, -1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) ** 2)),
|
||||
(0, np.inf, -1, np.inf, -1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) ** 2)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-1, inf] for all n.
|
||||
(-1, np.inf, -1, np.inf, -1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) ** 3)),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = [1, inf] and Dy = Dz = [-1, inf].
|
||||
(1, np.inf, -1, np.inf, -1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (((erf(1) + 1) ** 2) * erfc(1))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = Dy = [1, inf] and Dz = [-1, inf].
|
||||
(1, np.inf, 1, np.inf, -1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) * (erfc(1) ** 2))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = Dz = [1, inf] and Dy = [-1, inf].
|
||||
(1, np.inf, -1, np.inf, 1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) * (erfc(1) ** 2))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = [-1, inf] and Dy = Dz = [1, inf].
|
||||
(-1, np.inf, 1, np.inf, 1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * ((erf(1) + 1) * (erfc(1) ** 2))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = Dy = [-1, inf] and Dz = [1, inf].
|
||||
(-1, np.inf, -1, np.inf, 1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (((erf(1) + 1) ** 2) * erfc(1))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain Dx = Dz = [-1, inf] and Dy = [1, inf].
|
||||
(-1, np.inf, 1, np.inf, -1, np.inf,
|
||||
(np.pi ** (3 / 2)) / 8 * (((erf(1) + 1) ** 2) * erfc(1))),
|
||||
# Multiple integration of a function in n = 3 variables: f(x, y, z)
|
||||
# over domain D = [-inf, inf] for all n.
|
||||
(-np.inf, np.inf, -np.inf, np.inf, -np.inf, np.inf,
|
||||
np.pi ** (3 / 2)),
|
||||
],
|
||||
)
|
||||
def test_triple_integral_improper(
|
||||
self,
|
||||
x_lower,
|
||||
x_upper,
|
||||
y_lower,
|
||||
y_upper,
|
||||
z_lower,
|
||||
z_upper,
|
||||
expected
|
||||
):
|
||||
# The Gaussian Integral.
|
||||
def f(x, y, z):
|
||||
return np.exp(-x ** 2 - y ** 2 - z ** 2)
|
||||
|
||||
assert_quad(
|
||||
tplquad(f, x_lower, x_upper, y_lower, y_upper, z_lower, z_upper),
|
||||
expected,
|
||||
error_tolerance=6e-8
|
||||
)
|
||||
|
||||
def test_complex(self):
|
||||
def tfunc(x):
|
||||
return np.exp(1j*x)
|
||||
|
||||
assert np.allclose(
|
||||
quad(tfunc, 0, np.pi/2, complex_func=True)[0],
|
||||
1+1j)
|
||||
|
||||
# We consider a divergent case in order to force quadpack
|
||||
# to return an error message. The output is compared
|
||||
# against what is returned by explicit integration
|
||||
# of the parts.
|
||||
kwargs = {'a': 0, 'b': np.inf, 'full_output': True,
|
||||
'weight': 'cos', 'wvar': 1}
|
||||
res_c = quad(tfunc, complex_func=True, **kwargs)
|
||||
res_r = quad(lambda x: np.real(np.exp(1j*x)),
|
||||
complex_func=False,
|
||||
**kwargs)
|
||||
res_i = quad(lambda x: np.imag(np.exp(1j*x)),
|
||||
complex_func=False,
|
||||
**kwargs)
|
||||
|
||||
np.testing.assert_equal(res_c[0], res_r[0] + 1j*res_i[0])
|
||||
np.testing.assert_equal(res_c[1], res_r[1] + 1j*res_i[1])
|
||||
|
||||
assert len(res_c[2]['real']) == len(res_r[2:]) == 3
|
||||
assert res_c[2]['real'][2] == res_r[4]
|
||||
assert res_c[2]['real'][1] == res_r[3]
|
||||
assert res_c[2]['real'][0]['lst'] == res_r[2]['lst']
|
||||
|
||||
assert len(res_c[2]['imag']) == len(res_i[2:]) == 1
|
||||
assert res_c[2]['imag'][0]['lst'] == res_i[2]['lst']
|
||||
|
||||
|
||||
class TestNQuad:
|
||||
@pytest.mark.fail_slow(5)
|
||||
def test_fixed_limits(self):
|
||||
def func1(x0, x1, x2, x3):
|
||||
val = (x0**2 + x1*x2 - x3**3 + np.sin(x0) +
|
||||
(1 if (x0 - 0.2*x3 - 0.5 - 0.25*x1 > 0) else 0))
|
||||
return val
|
||||
|
||||
def opts_basic(*args):
|
||||
return {'points': [0.2*args[2] + 0.5 + 0.25*args[0]]}
|
||||
|
||||
res = nquad(func1, [[0, 1], [-1, 1], [.13, .8], [-.15, 1]],
|
||||
opts=[opts_basic, {}, {}, {}], full_output=True)
|
||||
assert_quad(res[:-1], 1.5267454070738635)
|
||||
assert_(res[-1]['neval'] > 0 and res[-1]['neval'] < 4e5)
|
||||
|
||||
@pytest.mark.fail_slow(5)
|
||||
def test_variable_limits(self):
|
||||
scale = .1
|
||||
|
||||
def func2(x0, x1, x2, x3, t0, t1):
|
||||
val = (x0*x1*x3**2 + np.sin(x2) + 1 +
|
||||
(1 if x0 + t1*x1 - t0 > 0 else 0))
|
||||
return val
|
||||
|
||||
def lim0(x1, x2, x3, t0, t1):
|
||||
return [scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) - 1,
|
||||
scale * (x1**2 + x2 + np.cos(x3)*t0*t1 + 1) + 1]
|
||||
|
||||
def lim1(x2, x3, t0, t1):
|
||||
return [scale * (t0*x2 + t1*x3) - 1,
|
||||
scale * (t0*x2 + t1*x3) + 1]
|
||||
|
||||
def lim2(x3, t0, t1):
|
||||
return [scale * (x3 + t0**2*t1**3) - 1,
|
||||
scale * (x3 + t0**2*t1**3) + 1]
|
||||
|
||||
def lim3(t0, t1):
|
||||
return [scale * (t0 + t1) - 1, scale * (t0 + t1) + 1]
|
||||
|
||||
def opts0(x1, x2, x3, t0, t1):
|
||||
return {'points': [t0 - t1*x1]}
|
||||
|
||||
def opts1(x2, x3, t0, t1):
|
||||
return {}
|
||||
|
||||
def opts2(x3, t0, t1):
|
||||
return {}
|
||||
|
||||
def opts3(t0, t1):
|
||||
return {}
|
||||
|
||||
res = nquad(func2, [lim0, lim1, lim2, lim3], args=(0, 0),
|
||||
opts=[opts0, opts1, opts2, opts3])
|
||||
assert_quad(res, 25.066666666666663)
|
||||
|
||||
def test_square_separate_ranges_and_opts(self):
|
||||
def f(y, x):
|
||||
return 1.0
|
||||
|
||||
assert_quad(nquad(f, [[-1, 1], [-1, 1]], opts=[{}, {}]), 4.0)
|
||||
|
||||
def test_square_aliased_ranges_and_opts(self):
|
||||
def f(y, x):
|
||||
return 1.0
|
||||
|
||||
r = [-1, 1]
|
||||
opt = {}
|
||||
assert_quad(nquad(f, [r, r], opts=[opt, opt]), 4.0)
|
||||
|
||||
def test_square_separate_fn_ranges_and_opts(self):
|
||||
def f(y, x):
|
||||
return 1.0
|
||||
|
||||
def fn_range0(*args):
|
||||
return (-1, 1)
|
||||
|
||||
def fn_range1(*args):
|
||||
return (-1, 1)
|
||||
|
||||
def fn_opt0(*args):
|
||||
return {}
|
||||
|
||||
def fn_opt1(*args):
|
||||
return {}
|
||||
|
||||
ranges = [fn_range0, fn_range1]
|
||||
opts = [fn_opt0, fn_opt1]
|
||||
assert_quad(nquad(f, ranges, opts=opts), 4.0)
|
||||
|
||||
def test_square_aliased_fn_ranges_and_opts(self):
|
||||
def f(y, x):
|
||||
return 1.0
|
||||
|
||||
def fn_range(*args):
|
||||
return (-1, 1)
|
||||
|
||||
def fn_opt(*args):
|
||||
return {}
|
||||
|
||||
ranges = [fn_range, fn_range]
|
||||
opts = [fn_opt, fn_opt]
|
||||
assert_quad(nquad(f, ranges, opts=opts), 4.0)
|
||||
|
||||
def test_matching_quad(self):
|
||||
def func(x):
|
||||
return x**2 + 1
|
||||
|
||||
res, reserr = quad(func, 0, 4)
|
||||
res2, reserr2 = nquad(func, ranges=[[0, 4]])
|
||||
assert_almost_equal(res, res2)
|
||||
assert_almost_equal(reserr, reserr2)
|
||||
|
||||
def test_matching_dblquad(self):
|
||||
def func2d(x0, x1):
|
||||
return x0**2 + x1**3 - x0 * x1 + 1
|
||||
|
||||
res, reserr = dblquad(func2d, -2, 2, lambda x: -3, lambda x: 3)
|
||||
res2, reserr2 = nquad(func2d, [[-3, 3], (-2, 2)])
|
||||
assert_almost_equal(res, res2)
|
||||
assert_almost_equal(reserr, reserr2)
|
||||
|
||||
def test_matching_tplquad(self):
|
||||
def func3d(x0, x1, x2, c0, c1):
|
||||
return x0**2 + c0 * x1**3 - x0 * x1 + 1 + c1 * np.sin(x2)
|
||||
|
||||
res = tplquad(func3d, -1, 2, lambda x: -2, lambda x: 2,
|
||||
lambda x, y: -np.pi, lambda x, y: np.pi,
|
||||
args=(2, 3))
|
||||
res2 = nquad(func3d, [[-np.pi, np.pi], [-2, 2], (-1, 2)], args=(2, 3))
|
||||
assert_almost_equal(res, res2)
|
||||
|
||||
def test_dict_as_opts(self):
|
||||
try:
|
||||
nquad(lambda x, y: x * y, [[0, 1], [0, 1]], opts={'epsrel': 0.0001})
|
||||
except TypeError:
|
||||
assert False
|
||||
|
||||
|
|
@ -1,730 +0,0 @@
|
|||
# mypy: disable-error-code="attr-defined"
|
||||
import pytest
|
||||
import numpy as np
|
||||
from numpy.testing import assert_equal, assert_almost_equal, assert_allclose
|
||||
from hypothesis import given
|
||||
import hypothesis.strategies as st
|
||||
import hypothesis.extra.numpy as hyp_num
|
||||
|
||||
from scipy.integrate import (romb, newton_cotes,
|
||||
cumulative_trapezoid, trapezoid,
|
||||
quad, simpson, fixed_quad,
|
||||
qmc_quad, cumulative_simpson)
|
||||
from scipy.integrate._quadrature import _cumulative_simpson_unequal_intervals
|
||||
|
||||
from scipy import stats, special, integrate
|
||||
from scipy.conftest import skip_xp_invalid_arg
|
||||
from scipy._lib._array_api_no_0d import xp_assert_close
|
||||
|
||||
skip_xp_backends = pytest.mark.skip_xp_backends
|
||||
|
||||
|
||||
class TestFixedQuad:
|
||||
def test_scalar(self):
|
||||
n = 4
|
||||
expected = 1/(2*n)
|
||||
got, _ = fixed_quad(lambda x: x**(2*n - 1), 0, 1, n=n)
|
||||
# quadrature exact for this input
|
||||
assert_allclose(got, expected, rtol=1e-12)
|
||||
|
||||
def test_vector(self):
|
||||
n = 4
|
||||
p = np.arange(1, 2*n)
|
||||
expected = 1/(p + 1)
|
||||
got, _ = fixed_quad(lambda x: x**p[:, None], 0, 1, n=n)
|
||||
assert_allclose(got, expected, rtol=1e-12)
|
||||
|
||||
|
||||
class TestQuadrature:
|
||||
def quad(self, x, a, b, args):
|
||||
raise NotImplementedError
|
||||
|
||||
def test_romb(self):
|
||||
assert_equal(romb(np.arange(17)), 128)
|
||||
|
||||
def test_romb_gh_3731(self):
|
||||
# Check that romb makes maximal use of data points
|
||||
x = np.arange(2**4+1)
|
||||
y = np.cos(0.2*x)
|
||||
val = romb(y)
|
||||
val2, err = quad(lambda x: np.cos(0.2*x), x.min(), x.max())
|
||||
assert_allclose(val, val2, rtol=1e-8, atol=0)
|
||||
|
||||
def test_newton_cotes(self):
|
||||
"""Test the first few degrees, for evenly spaced points."""
|
||||
n = 1
|
||||
wts, errcoff = newton_cotes(n, 1)
|
||||
assert_equal(wts, n*np.array([0.5, 0.5]))
|
||||
assert_almost_equal(errcoff, -n**3/12.0)
|
||||
|
||||
n = 2
|
||||
wts, errcoff = newton_cotes(n, 1)
|
||||
assert_almost_equal(wts, n*np.array([1.0, 4.0, 1.0])/6.0)
|
||||
assert_almost_equal(errcoff, -n**5/2880.0)
|
||||
|
||||
n = 3
|
||||
wts, errcoff = newton_cotes(n, 1)
|
||||
assert_almost_equal(wts, n*np.array([1.0, 3.0, 3.0, 1.0])/8.0)
|
||||
assert_almost_equal(errcoff, -n**5/6480.0)
|
||||
|
||||
n = 4
|
||||
wts, errcoff = newton_cotes(n, 1)
|
||||
assert_almost_equal(wts, n*np.array([7.0, 32.0, 12.0, 32.0, 7.0])/90.0)
|
||||
assert_almost_equal(errcoff, -n**7/1935360.0)
|
||||
|
||||
def test_newton_cotes2(self):
|
||||
"""Test newton_cotes with points that are not evenly spaced."""
|
||||
|
||||
x = np.array([0.0, 1.5, 2.0])
|
||||
y = x**2
|
||||
wts, errcoff = newton_cotes(x)
|
||||
exact_integral = 8.0/3
|
||||
numeric_integral = np.dot(wts, y)
|
||||
assert_almost_equal(numeric_integral, exact_integral)
|
||||
|
||||
x = np.array([0.0, 1.4, 2.1, 3.0])
|
||||
y = x**2
|
||||
wts, errcoff = newton_cotes(x)
|
||||
exact_integral = 9.0
|
||||
numeric_integral = np.dot(wts, y)
|
||||
assert_almost_equal(numeric_integral, exact_integral)
|
||||
|
||||
def test_simpson(self):
|
||||
y = np.arange(17)
|
||||
assert_equal(simpson(y), 128)
|
||||
assert_equal(simpson(y, dx=0.5), 64)
|
||||
assert_equal(simpson(y, x=np.linspace(0, 4, 17)), 32)
|
||||
|
||||
# integral should be exactly 21
|
||||
x = np.linspace(1, 4, 4)
|
||||
def f(x):
|
||||
return x**2
|
||||
|
||||
assert_allclose(simpson(f(x), x=x), 21.0)
|
||||
|
||||
# integral should be exactly 114
|
||||
x = np.linspace(1, 7, 4)
|
||||
assert_allclose(simpson(f(x), dx=2.0), 114)
|
||||
|
||||
# test multi-axis behaviour
|
||||
a = np.arange(16).reshape(4, 4)
|
||||
x = np.arange(64.).reshape(4, 4, 4)
|
||||
y = f(x)
|
||||
for i in range(3):
|
||||
r = simpson(y, x=x, axis=i)
|
||||
it = np.nditer(a, flags=['multi_index'])
|
||||
for _ in it:
|
||||
idx = list(it.multi_index)
|
||||
idx.insert(i, slice(None))
|
||||
integral = x[tuple(idx)][-1]**3 / 3 - x[tuple(idx)][0]**3 / 3
|
||||
assert_allclose(r[it.multi_index], integral)
|
||||
|
||||
# test when integration axis only has two points
|
||||
x = np.arange(16).reshape(8, 2)
|
||||
y = f(x)
|
||||
r = simpson(y, x=x, axis=-1)
|
||||
|
||||
integral = 0.5 * (y[:, 1] + y[:, 0]) * (x[:, 1] - x[:, 0])
|
||||
assert_allclose(r, integral)
|
||||
|
||||
# odd points, test multi-axis behaviour
|
||||
a = np.arange(25).reshape(5, 5)
|
||||
x = np.arange(125).reshape(5, 5, 5)
|
||||
y = f(x)
|
||||
for i in range(3):
|
||||
r = simpson(y, x=x, axis=i)
|
||||
it = np.nditer(a, flags=['multi_index'])
|
||||
for _ in it:
|
||||
idx = list(it.multi_index)
|
||||
idx.insert(i, slice(None))
|
||||
integral = x[tuple(idx)][-1]**3 / 3 - x[tuple(idx)][0]**3 / 3
|
||||
assert_allclose(r[it.multi_index], integral)
|
||||
|
||||
# Tests for checking base case
|
||||
x = np.array([3])
|
||||
y = np.power(x, 2)
|
||||
assert_allclose(simpson(y, x=x, axis=0), 0.0)
|
||||
assert_allclose(simpson(y, x=x, axis=-1), 0.0)
|
||||
|
||||
x = np.array([3, 3, 3, 3])
|
||||
y = np.power(x, 2)
|
||||
assert_allclose(simpson(y, x=x, axis=0), 0.0)
|
||||
assert_allclose(simpson(y, x=x, axis=-1), 0.0)
|
||||
|
||||
x = np.array([[1, 2, 4, 8], [1, 2, 4, 8], [1, 2, 4, 8]])
|
||||
y = np.power(x, 2)
|
||||
zero_axis = [0.0, 0.0, 0.0, 0.0]
|
||||
default_axis = [170 + 1/3] * 3 # 8**3 / 3 - 1/3
|
||||
assert_allclose(simpson(y, x=x, axis=0), zero_axis)
|
||||
# the following should be exact
|
||||
assert_allclose(simpson(y, x=x, axis=-1), default_axis)
|
||||
|
||||
x = np.array([[1, 2, 4, 8], [1, 2, 4, 8], [1, 8, 16, 32]])
|
||||
y = np.power(x, 2)
|
||||
zero_axis = [0.0, 136.0, 1088.0, 8704.0]
|
||||
default_axis = [170 + 1/3, 170 + 1/3, 32**3 / 3 - 1/3]
|
||||
assert_allclose(simpson(y, x=x, axis=0), zero_axis)
|
||||
assert_allclose(simpson(y, x=x, axis=-1), default_axis)
|
||||
|
||||
|
||||
@pytest.mark.parametrize('droplast', [False, True])
|
||||
def test_simpson_2d_integer_no_x(self, droplast):
|
||||
# The inputs are 2d integer arrays. The results should be
|
||||
# identical to the results when the inputs are floating point.
|
||||
y = np.array([[2, 2, 4, 4, 8, 8, -4, 5],
|
||||
[4, 4, 2, -4, 10, 22, -2, 10]])
|
||||
if droplast:
|
||||
y = y[:, :-1]
|
||||
result = simpson(y, axis=-1)
|
||||
expected = simpson(np.array(y, dtype=np.float64), axis=-1)
|
||||
assert_equal(result, expected)
|
||||
|
||||
|
||||
class TestCumulative_trapezoid:
|
||||
def test_1d(self):
|
||||
x = np.linspace(-2, 2, num=5)
|
||||
y = x
|
||||
y_int = cumulative_trapezoid(y, x, initial=0)
|
||||
y_expected = [0., -1.5, -2., -1.5, 0.]
|
||||
assert_allclose(y_int, y_expected)
|
||||
|
||||
y_int = cumulative_trapezoid(y, x, initial=None)
|
||||
assert_allclose(y_int, y_expected[1:])
|
||||
|
||||
def test_y_nd_x_nd(self):
|
||||
x = np.arange(3 * 2 * 4).reshape(3, 2, 4)
|
||||
y = x
|
||||
y_int = cumulative_trapezoid(y, x, initial=0)
|
||||
y_expected = np.array([[[0., 0.5, 2., 4.5],
|
||||
[0., 4.5, 10., 16.5]],
|
||||
[[0., 8.5, 18., 28.5],
|
||||
[0., 12.5, 26., 40.5]],
|
||||
[[0., 16.5, 34., 52.5],
|
||||
[0., 20.5, 42., 64.5]]])
|
||||
|
||||
assert_allclose(y_int, y_expected)
|
||||
|
||||
# Try with all axes
|
||||
shapes = [(2, 2, 4), (3, 1, 4), (3, 2, 3)]
|
||||
for axis, shape in zip([0, 1, 2], shapes):
|
||||
y_int = cumulative_trapezoid(y, x, initial=0, axis=axis)
|
||||
assert_equal(y_int.shape, (3, 2, 4))
|
||||
y_int = cumulative_trapezoid(y, x, initial=None, axis=axis)
|
||||
assert_equal(y_int.shape, shape)
|
||||
|
||||
def test_y_nd_x_1d(self):
|
||||
y = np.arange(3 * 2 * 4).reshape(3, 2, 4)
|
||||
x = np.arange(4)**2
|
||||
# Try with all axes
|
||||
ys_expected = (
|
||||
np.array([[[4., 5., 6., 7.],
|
||||
[8., 9., 10., 11.]],
|
||||
[[40., 44., 48., 52.],
|
||||
[56., 60., 64., 68.]]]),
|
||||
np.array([[[2., 3., 4., 5.]],
|
||||
[[10., 11., 12., 13.]],
|
||||
[[18., 19., 20., 21.]]]),
|
||||
np.array([[[0.5, 5., 17.5],
|
||||
[4.5, 21., 53.5]],
|
||||
[[8.5, 37., 89.5],
|
||||
[12.5, 53., 125.5]],
|
||||
[[16.5, 69., 161.5],
|
||||
[20.5, 85., 197.5]]]))
|
||||
|
||||
for axis, y_expected in zip([0, 1, 2], ys_expected):
|
||||
y_int = cumulative_trapezoid(y, x=x[:y.shape[axis]], axis=axis,
|
||||
initial=None)
|
||||
assert_allclose(y_int, y_expected)
|
||||
|
||||
def test_x_none(self):
|
||||
y = np.linspace(-2, 2, num=5)
|
||||
|
||||
y_int = cumulative_trapezoid(y)
|
||||
y_expected = [-1.5, -2., -1.5, 0.]
|
||||
assert_allclose(y_int, y_expected)
|
||||
|
||||
y_int = cumulative_trapezoid(y, initial=0)
|
||||
y_expected = [0, -1.5, -2., -1.5, 0.]
|
||||
assert_allclose(y_int, y_expected)
|
||||
|
||||
y_int = cumulative_trapezoid(y, dx=3)
|
||||
y_expected = [-4.5, -6., -4.5, 0.]
|
||||
assert_allclose(y_int, y_expected)
|
||||
|
||||
y_int = cumulative_trapezoid(y, dx=3, initial=0)
|
||||
y_expected = [0, -4.5, -6., -4.5, 0.]
|
||||
assert_allclose(y_int, y_expected)
|
||||
|
||||
@pytest.mark.parametrize(
|
||||
"initial", [1, 0.5]
|
||||
)
|
||||
def test_initial_error(self, initial):
|
||||
"""If initial is not None or 0, a ValueError is raised."""
|
||||
y = np.linspace(0, 10, num=10)
|
||||
with pytest.raises(ValueError, match="`initial`"):
|
||||
cumulative_trapezoid(y, initial=initial)
|
||||
|
||||
def test_zero_len_y(self):
|
||||
with pytest.raises(ValueError, match="At least one point is required"):
|
||||
cumulative_trapezoid(y=[])
|
||||
|
||||
|
||||
class TestTrapezoid:
|
||||
def test_simple(self, xp):
|
||||
x = xp.arange(-10, 10, .1)
|
||||
r = trapezoid(xp.exp(-.5 * x ** 2) / xp.sqrt(2 * xp.asarray(xp.pi)), dx=0.1)
|
||||
# check integral of normal equals 1
|
||||
xp_assert_close(r, xp.asarray(1.0))
|
||||
|
||||
@skip_xp_backends('jax.numpy',
|
||||
reason="JAX arrays do not support item assignment")
|
||||
def test_ndim(self, xp):
|
||||
x = xp.linspace(0, 1, 3)
|
||||
y = xp.linspace(0, 2, 8)
|
||||
z = xp.linspace(0, 3, 13)
|
||||
|
||||
wx = xp.ones_like(x) * (x[1] - x[0])
|
||||
wx[0] /= 2
|
||||
wx[-1] /= 2
|
||||
wy = xp.ones_like(y) * (y[1] - y[0])
|
||||
wy[0] /= 2
|
||||
wy[-1] /= 2
|
||||
wz = xp.ones_like(z) * (z[1] - z[0])
|
||||
wz[0] /= 2
|
||||
wz[-1] /= 2
|
||||
|
||||
q = x[:, None, None] + y[None,:, None] + z[None, None,:]
|
||||
|
||||
qx = xp.sum(q * wx[:, None, None], axis=0)
|
||||
qy = xp.sum(q * wy[None, :, None], axis=1)
|
||||
qz = xp.sum(q * wz[None, None, :], axis=2)
|
||||
|
||||
# n-d `x`
|
||||
r = trapezoid(q, x=x[:, None, None], axis=0)
|
||||
xp_assert_close(r, qx)
|
||||
r = trapezoid(q, x=y[None,:, None], axis=1)
|
||||
xp_assert_close(r, qy)
|
||||
r = trapezoid(q, x=z[None, None,:], axis=2)
|
||||
xp_assert_close(r, qz)
|
||||
|
||||
# 1-d `x`
|
||||
r = trapezoid(q, x=x, axis=0)
|
||||
xp_assert_close(r, qx)
|
||||
r = trapezoid(q, x=y, axis=1)
|
||||
xp_assert_close(r, qy)
|
||||
r = trapezoid(q, x=z, axis=2)
|
||||
xp_assert_close(r, qz)
|
||||
|
||||
@skip_xp_backends('jax.numpy',
|
||||
reason="JAX arrays do not support item assignment")
|
||||
def test_gh21908(self, xp):
|
||||
# extended testing for n-dim arrays
|
||||
x = xp.reshape(xp.linspace(0, 29, 30), (3, 10))
|
||||
y = xp.reshape(xp.linspace(0, 29, 30), (3, 10))
|
||||
|
||||
out0 = xp.linspace(200, 380, 10)
|
||||
xp_assert_close(trapezoid(y, x=x, axis=0), out0)
|
||||
xp_assert_close(trapezoid(y, x=xp.asarray([0, 10., 20.]), axis=0), out0)
|
||||
# x needs to be broadcastable against y
|
||||
xp_assert_close(
|
||||
trapezoid(y, x=xp.asarray([0, 10., 20.])[:, None], axis=0),
|
||||
out0
|
||||
)
|
||||
with pytest.raises(Exception):
|
||||
# x is not broadcastable against y
|
||||
trapezoid(y, x=xp.asarray([0, 10., 20.])[None, :], axis=0)
|
||||
|
||||
out1 = xp.asarray([ 40.5, 130.5, 220.5])
|
||||
xp_assert_close(trapezoid(y, x=x, axis=1), out1)
|
||||
xp_assert_close(
|
||||
trapezoid(y, x=xp.linspace(0, 9, 10), axis=1),
|
||||
out1
|
||||
)
|
||||
|
||||
@skip_xp_invalid_arg
|
||||
def test_masked(self, xp):
|
||||
# Testing that masked arrays behave as if the function is 0 where
|
||||
# masked
|
||||
x = np.arange(5)
|
||||
y = x * x
|
||||
mask = x == 2
|
||||
ym = np.ma.array(y, mask=mask)
|
||||
r = 13.0 # sum(0.5 * (0 + 1) * 1.0 + 0.5 * (9 + 16))
|
||||
assert_allclose(trapezoid(ym, x), r)
|
||||
|
||||
xm = np.ma.array(x, mask=mask)
|
||||
assert_allclose(trapezoid(ym, xm), r)
|
||||
|
||||
xm = np.ma.array(x, mask=mask)
|
||||
assert_allclose(trapezoid(y, xm), r)
|
||||
|
||||
@skip_xp_backends(np_only=True,
|
||||
reason='array-likes only supported for NumPy backend')
|
||||
def test_array_like(self, xp):
|
||||
x = list(range(5))
|
||||
y = [t * t for t in x]
|
||||
xarr = xp.asarray(x, dtype=xp.float64)
|
||||
yarr = xp.asarray(y, dtype=xp.float64)
|
||||
res = trapezoid(y, x)
|
||||
resarr = trapezoid(yarr, xarr)
|
||||
xp_assert_close(res, resarr)
|
||||
|
||||
|
||||
class TestQMCQuad:
|
||||
@pytest.mark.thread_unsafe
|
||||
def test_input_validation(self):
|
||||
message = "`func` must be callable."
|
||||
with pytest.raises(TypeError, match=message):
|
||||
qmc_quad("a duck", [0, 0], [1, 1])
|
||||
|
||||
message = "`func` must evaluate the integrand at points..."
|
||||
with pytest.raises(ValueError, match=message):
|
||||
qmc_quad(lambda: 1, [0, 0], [1, 1])
|
||||
|
||||
def func(x):
|
||||
assert x.ndim == 1
|
||||
return np.sum(x)
|
||||
message = "Exception encountered when attempting vectorized call..."
|
||||
with pytest.warns(UserWarning, match=message):
|
||||
qmc_quad(func, [0, 0], [1, 1])
|
||||
|
||||
message = "`n_points` must be an integer."
|
||||
with pytest.raises(TypeError, match=message):
|
||||
qmc_quad(lambda x: 1, [0, 0], [1, 1], n_points=1024.5)
|
||||
|
||||
message = "`n_estimates` must be an integer."
|
||||
with pytest.raises(TypeError, match=message):
|
||||
qmc_quad(lambda x: 1, [0, 0], [1, 1], n_estimates=8.5)
|
||||
|
||||
message = "`qrng` must be an instance of scipy.stats.qmc.QMCEngine."
|
||||
with pytest.raises(TypeError, match=message):
|
||||
qmc_quad(lambda x: 1, [0, 0], [1, 1], qrng="a duck")
|
||||
|
||||
message = "`qrng` must be initialized with dimensionality equal to "
|
||||
with pytest.raises(ValueError, match=message):
|
||||
qmc_quad(lambda x: 1, [0, 0], [1, 1], qrng=stats.qmc.Sobol(1))
|
||||
|
||||
message = r"`log` must be boolean \(`True` or `False`\)."
|
||||
with pytest.raises(TypeError, match=message):
|
||||
qmc_quad(lambda x: 1, [0, 0], [1, 1], log=10)
|
||||
|
||||
def basic_test(self, n_points=2**8, n_estimates=8, signs=None):
|
||||
if signs is None:
|
||||
signs = np.ones(2)
|
||||
ndim = 2
|
||||
mean = np.zeros(ndim)
|
||||
cov = np.eye(ndim)
|
||||
|
||||
def func(x):
|
||||
return stats.multivariate_normal.pdf(x.T, mean, cov)
|
||||
|
||||
rng = np.random.default_rng(2879434385674690281)
|
||||
qrng = stats.qmc.Sobol(ndim, seed=rng)
|
||||
a = np.zeros(ndim)
|
||||
b = np.ones(ndim) * signs
|
||||
res = qmc_quad(func, a, b, n_points=n_points,
|
||||
n_estimates=n_estimates, qrng=qrng)
|
||||
ref = stats.multivariate_normal.cdf(b, mean, cov, lower_limit=a)
|
||||
atol = special.stdtrit(n_estimates-1, 0.995) * res.standard_error # 99% CI
|
||||
assert_allclose(res.integral, ref, atol=atol)
|
||||
assert np.prod(signs)*res.integral > 0
|
||||
|
||||
rng = np.random.default_rng(2879434385674690281)
|
||||
qrng = stats.qmc.Sobol(ndim, seed=rng)
|
||||
logres = qmc_quad(lambda *args: np.log(func(*args)), a, b,
|
||||
n_points=n_points, n_estimates=n_estimates,
|
||||
log=True, qrng=qrng)
|
||||
assert_allclose(np.exp(logres.integral), res.integral, rtol=1e-14)
|
||||
assert np.imag(logres.integral) == (np.pi if np.prod(signs) < 0 else 0)
|
||||
assert_allclose(np.exp(logres.standard_error),
|
||||
res.standard_error, rtol=1e-14, atol=1e-16)
|
||||
|
||||
@pytest.mark.parametrize("n_points", [2**8, 2**12])
|
||||
@pytest.mark.parametrize("n_estimates", [8, 16])
|
||||
def test_basic(self, n_points, n_estimates):
|
||||
self.basic_test(n_points, n_estimates)
|
||||
|
||||
@pytest.mark.parametrize("signs", [[1, 1], [-1, -1], [-1, 1], [1, -1]])
|
||||
def test_sign(self, signs):
|
||||
self.basic_test(signs=signs)
|
||||
|
||||
@pytest.mark.thread_unsafe
|
||||
@pytest.mark.parametrize("log", [False, True])
|
||||
def test_zero(self, log):
|
||||
message = "A lower limit was equal to an upper limit, so"
|
||||
with pytest.warns(UserWarning, match=message):
|
||||
res = qmc_quad(lambda x: 1, [0, 0], [0, 1], log=log)
|
||||
assert res.integral == (-np.inf if log else 0)
|
||||
assert res.standard_error == 0
|
||||
|
||||
def test_flexible_input(self):
|
||||
# check that qrng is not required
|
||||
# also checks that for 1d problems, a and b can be scalars
|
||||
def func(x):
|
||||
return stats.norm.pdf(x, scale=2)
|
||||
|
||||
res = qmc_quad(func, 0, 1)
|
||||
ref = stats.norm.cdf(1, scale=2) - stats.norm.cdf(0, scale=2)
|
||||
assert_allclose(res.integral, ref, 1e-2)
|
||||
|
||||
|
||||
def cumulative_simpson_nd_reference(y, *, x=None, dx=None, initial=None, axis=-1):
|
||||
# Use cumulative_trapezoid if length of y < 3
|
||||
if y.shape[axis] < 3:
|
||||
if initial is None:
|
||||
return cumulative_trapezoid(y, x=x, dx=dx, axis=axis, initial=None)
|
||||
else:
|
||||
return initial + cumulative_trapezoid(y, x=x, dx=dx, axis=axis, initial=0)
|
||||
|
||||
# Ensure that working axis is last axis
|
||||
y = np.moveaxis(y, axis, -1)
|
||||
x = np.moveaxis(x, axis, -1) if np.ndim(x) > 1 else x
|
||||
dx = np.moveaxis(dx, axis, -1) if np.ndim(dx) > 1 else dx
|
||||
initial = np.moveaxis(initial, axis, -1) if np.ndim(initial) > 1 else initial
|
||||
|
||||
# If `x` is not present, create it from `dx`
|
||||
n = y.shape[-1]
|
||||
x = dx * np.arange(n) if dx is not None else x
|
||||
# Similarly, if `initial` is not present, set it to 0
|
||||
initial_was_none = initial is None
|
||||
initial = 0 if initial_was_none else initial
|
||||
|
||||
# `np.apply_along_axis` accepts only one array, so concatenate arguments
|
||||
x = np.broadcast_to(x, y.shape)
|
||||
initial = np.broadcast_to(initial, y.shape[:-1] + (1,))
|
||||
z = np.concatenate((y, x, initial), axis=-1)
|
||||
|
||||
# Use `np.apply_along_axis` to compute result
|
||||
def f(z):
|
||||
return cumulative_simpson(z[:n], x=z[n:2*n], initial=z[2*n:])
|
||||
res = np.apply_along_axis(f, -1, z)
|
||||
|
||||
# Remove `initial` and undo axis move as needed
|
||||
res = res[..., 1:] if initial_was_none else res
|
||||
res = np.moveaxis(res, -1, axis)
|
||||
return res
|
||||
|
||||
|
||||
class TestCumulativeSimpson:
|
||||
x0 = np.arange(4)
|
||||
y0 = x0**2
|
||||
|
||||
@pytest.mark.parametrize('use_dx', (False, True))
|
||||
@pytest.mark.parametrize('use_initial', (False, True))
|
||||
def test_1d(self, use_dx, use_initial):
|
||||
# Test for exact agreement with polynomial of highest
|
||||
# possible order (3 if `dx` is constant, 2 otherwise).
|
||||
rng = np.random.default_rng(82456839535679456794)
|
||||
n = 10
|
||||
|
||||
# Generate random polynomials and ground truth
|
||||
# integral of appropriate order
|
||||
order = 3 if use_dx else 2
|
||||
dx = rng.random()
|
||||
x = (np.sort(rng.random(n)) if order == 2
|
||||
else np.arange(n)*dx + rng.random())
|
||||
i = np.arange(order + 1)[:, np.newaxis]
|
||||
c = rng.random(order + 1)[:, np.newaxis]
|
||||
y = np.sum(c*x**i, axis=0)
|
||||
Y = np.sum(c*x**(i + 1)/(i + 1), axis=0)
|
||||
ref = Y if use_initial else (Y-Y[0])[1:]
|
||||
|
||||
# Integrate with `cumulative_simpson`
|
||||
initial = Y[0] if use_initial else None
|
||||
kwarg = {'dx': dx} if use_dx else {'x': x}
|
||||
res = cumulative_simpson(y, **kwarg, initial=initial)
|
||||
|
||||
# Compare result against reference
|
||||
if not use_dx:
|
||||
assert_allclose(res, ref, rtol=2e-15)
|
||||
else:
|
||||
i0 = 0 if use_initial else 1
|
||||
# all terms are "close"
|
||||
assert_allclose(res, ref, rtol=0.0025)
|
||||
# only even-interval terms are "exact"
|
||||
assert_allclose(res[i0::2], ref[i0::2], rtol=2e-15)
|
||||
|
||||
@pytest.mark.parametrize('axis', np.arange(-3, 3))
|
||||
@pytest.mark.parametrize('x_ndim', (1, 3))
|
||||
@pytest.mark.parametrize('x_len', (1, 2, 7))
|
||||
@pytest.mark.parametrize('i_ndim', (None, 0, 3,))
|
||||
@pytest.mark.parametrize('dx', (None, True))
|
||||
def test_nd(self, axis, x_ndim, x_len, i_ndim, dx):
|
||||
# Test behavior of `cumulative_simpson` with N-D `y`
|
||||
rng = np.random.default_rng(82456839535679456794)
|
||||
|
||||
# determine shapes
|
||||
shape = [5, 6, x_len]
|
||||
shape[axis], shape[-1] = shape[-1], shape[axis]
|
||||
shape_len_1 = shape.copy()
|
||||
shape_len_1[axis] = 1
|
||||
i_shape = shape_len_1 if i_ndim == 3 else ()
|
||||
|
||||
# initialize arguments
|
||||
y = rng.random(size=shape)
|
||||
x, dx = None, None
|
||||
if dx:
|
||||
dx = rng.random(size=shape_len_1) if x_ndim > 1 else rng.random()
|
||||
else:
|
||||
x = (np.sort(rng.random(size=shape), axis=axis) if x_ndim > 1
|
||||
else np.sort(rng.random(size=shape[axis])))
|
||||
initial = None if i_ndim is None else rng.random(size=i_shape)
|
||||
|
||||
# compare results
|
||||
res = cumulative_simpson(y, x=x, dx=dx, initial=initial, axis=axis)
|
||||
ref = cumulative_simpson_nd_reference(y, x=x, dx=dx, initial=initial, axis=axis)
|
||||
np.testing.assert_allclose(res, ref, rtol=1e-15)
|
||||
|
||||
@pytest.mark.parametrize(('message', 'kwarg_update'), [
|
||||
("x must be strictly increasing", dict(x=[2, 2, 3, 4])),
|
||||
("x must be strictly increasing", dict(x=[x0, [2, 2, 4, 8]], y=[y0, y0])),
|
||||
("x must be strictly increasing", dict(x=[x0, x0, x0], y=[y0, y0, y0], axis=0)),
|
||||
("At least one point is required", dict(x=[], y=[])),
|
||||
("`axis=4` is not valid for `y` with `y.ndim=1`", dict(axis=4)),
|
||||
("shape of `x` must be the same as `y` or 1-D", dict(x=np.arange(5))),
|
||||
("`initial` must either be a scalar or...", dict(initial=np.arange(5))),
|
||||
("`dx` must either be a scalar or...", dict(x=None, dx=np.arange(5))),
|
||||
])
|
||||
def test_simpson_exceptions(self, message, kwarg_update):
|
||||
kwargs0 = dict(y=self.y0, x=self.x0, dx=None, initial=None, axis=-1)
|
||||
with pytest.raises(ValueError, match=message):
|
||||
cumulative_simpson(**dict(kwargs0, **kwarg_update))
|
||||
|
||||
def test_special_cases(self):
|
||||
# Test special cases not checked elsewhere
|
||||
rng = np.random.default_rng(82456839535679456794)
|
||||
y = rng.random(size=10)
|
||||
res = cumulative_simpson(y, dx=0)
|
||||
assert_equal(res, 0)
|
||||
|
||||
# Should add tests of:
|
||||
# - all elements of `x` identical
|
||||
# These should work as they do for `simpson`
|
||||
|
||||
def _get_theoretical_diff_between_simps_and_cum_simps(self, y, x):
|
||||
"""`cumulative_simpson` and `simpson` can be tested against other to verify
|
||||
they give consistent results. `simpson` will iteratively be called with
|
||||
successively higher upper limits of integration. This function calculates
|
||||
the theoretical correction required to `simpson` at even intervals to match
|
||||
with `cumulative_simpson`.
|
||||
"""
|
||||
d = np.diff(x, axis=-1)
|
||||
sub_integrals_h1 = _cumulative_simpson_unequal_intervals(y, d)
|
||||
sub_integrals_h2 = _cumulative_simpson_unequal_intervals(
|
||||
y[..., ::-1], d[..., ::-1]
|
||||
)[..., ::-1]
|
||||
|
||||
# Concatenate to build difference array
|
||||
zeros_shape = (*y.shape[:-1], 1)
|
||||
theoretical_difference = np.concatenate(
|
||||
[
|
||||
np.zeros(zeros_shape),
|
||||
(sub_integrals_h1[..., 1:] - sub_integrals_h2[..., :-1]),
|
||||
np.zeros(zeros_shape),
|
||||
],
|
||||
axis=-1,
|
||||
)
|
||||
# Differences only expected at even intervals. Odd intervals will
|
||||
# match exactly so there is no correction
|
||||
theoretical_difference[..., 1::2] = 0.0
|
||||
# Note: the first interval will not match from this correction as
|
||||
# `simpson` uses the trapezoidal rule
|
||||
return theoretical_difference
|
||||
|
||||
@pytest.mark.fail_slow(10)
|
||||
@pytest.mark.thread_unsafe
|
||||
@pytest.mark.slow
|
||||
@given(
|
||||
y=hyp_num.arrays(
|
||||
np.float64,
|
||||
hyp_num.array_shapes(max_dims=4, min_side=3, max_side=10),
|
||||
elements=st.floats(-10, 10, allow_nan=False).filter(lambda x: abs(x) > 1e-7)
|
||||
)
|
||||
)
|
||||
def test_cumulative_simpson_against_simpson_with_default_dx(
|
||||
self, y
|
||||
):
|
||||
"""Theoretically, the output of `cumulative_simpson` will be identical
|
||||
to `simpson` at all even indices and in the last index. The first index
|
||||
will not match as `simpson` uses the trapezoidal rule when there are only two
|
||||
data points. Odd indices after the first index are shown to match with
|
||||
a mathematically-derived correction."""
|
||||
def simpson_reference(y):
|
||||
return np.stack(
|
||||
[simpson(y[..., :i], dx=1.0) for i in range(2, y.shape[-1]+1)], axis=-1,
|
||||
)
|
||||
|
||||
res = cumulative_simpson(y, dx=1.0)
|
||||
ref = simpson_reference(y)
|
||||
theoretical_difference = self._get_theoretical_diff_between_simps_and_cum_simps(
|
||||
y, x=np.arange(y.shape[-1])
|
||||
)
|
||||
np.testing.assert_allclose(
|
||||
res[..., 1:], ref[..., 1:] + theoretical_difference[..., 1:], atol=1e-16
|
||||
)
|
||||
|
||||
@pytest.mark.fail_slow(10)
|
||||
@pytest.mark.thread_unsafe
|
||||
@pytest.mark.slow
|
||||
@given(
|
||||
y=hyp_num.arrays(
|
||||
np.float64,
|
||||
hyp_num.array_shapes(max_dims=4, min_side=3, max_side=10),
|
||||
elements=st.floats(-10, 10, allow_nan=False).filter(lambda x: abs(x) > 1e-7)
|
||||
)
|
||||
)
|
||||
def test_cumulative_simpson_against_simpson(
|
||||
self, y
|
||||
):
|
||||
"""Theoretically, the output of `cumulative_simpson` will be identical
|
||||
to `simpson` at all even indices and in the last index. The first index
|
||||
will not match as `simpson` uses the trapezoidal rule when there are only two
|
||||
data points. Odd indices after the first index are shown to match with
|
||||
a mathematically-derived correction."""
|
||||
interval = 10/(y.shape[-1] - 1)
|
||||
x = np.linspace(0, 10, num=y.shape[-1])
|
||||
x[1:] = x[1:] + 0.2*interval*np.random.uniform(-1, 1, len(x) - 1)
|
||||
|
||||
def simpson_reference(y, x):
|
||||
return np.stack(
|
||||
[simpson(y[..., :i], x=x[..., :i]) for i in range(2, y.shape[-1]+1)],
|
||||
axis=-1,
|
||||
)
|
||||
|
||||
res = cumulative_simpson(y, x=x)
|
||||
ref = simpson_reference(y, x)
|
||||
theoretical_difference = self._get_theoretical_diff_between_simps_and_cum_simps(
|
||||
y, x
|
||||
)
|
||||
np.testing.assert_allclose(
|
||||
res[..., 1:], ref[..., 1:] + theoretical_difference[..., 1:]
|
||||
)
|
||||
|
||||
class TestLebedev:
|
||||
def test_input_validation(self):
|
||||
# only certain rules are available
|
||||
message = "Order n=-1 not available..."
|
||||
with pytest.raises(NotImplementedError, match=message):
|
||||
integrate.lebedev_rule(-1)
|
||||
|
||||
def test_quadrature(self):
|
||||
# Test points/weights to integrate an example function
|
||||
|
||||
def f(x):
|
||||
return np.exp(x[0])
|
||||
|
||||
x, w = integrate.lebedev_rule(15)
|
||||
res = w @ f(x)
|
||||
ref = 14.7680137457653 # lebedev_rule reference [3]
|
||||
assert_allclose(res, ref, rtol=1e-14)
|
||||
assert_allclose(np.sum(w), 4 * np.pi)
|
||||
|
||||
@pytest.mark.parametrize('order', list(range(3, 32, 2)) + list(range(35, 132, 6)))
|
||||
def test_properties(self, order):
|
||||
x, w = integrate.lebedev_rule(order)
|
||||
# dispersion should be maximal; no clear spherical mean
|
||||
with np.errstate(divide='ignore', invalid='ignore'):
|
||||
res = stats.directional_stats(x.T, axis=0)
|
||||
assert_allclose(res.mean_resultant_length, 0, atol=1e-15)
|
||||
# weights should sum to 4*pi (surface area of unit sphere)
|
||||
assert_allclose(np.sum(w), 4*np.pi)
|
||||
File diff suppressed because it is too large
Load diff
Loading…
Add table
Add a link
Reference in a new issue