Skip to content

Commit

Permalink
Dense output for yp.
Browse files Browse the repository at this point in the history
  • Loading branch information
JonasBreuling committed Aug 2, 2024
1 parent 0003086 commit fa0c9ee
Show file tree
Hide file tree
Showing 8 changed files with 373 additions and 88 deletions.
18 changes: 8 additions & 10 deletions examples/particle_on_circular_track.py
Original file line number Diff line number Diff line change
Expand Up @@ -55,8 +55,9 @@ def sol_true(t):
if __name__ == "__main__":
# time span
t0 = 1
t1 = t0 + 10
t1 = t0 + 6
t_span = (t0, t1)
t_eval = np.linspace(t0, t1, num=int(1e3))

# method = "BDF"
method = "Radau"
Expand All @@ -79,13 +80,12 @@ def sol_true(t):
# dae solution
##############
start = time.time()
sol = solve_dae(F, t_span, y0, yp0, atol=atol, rtol=rtol, method=method)
sol = solve_dae(F, t_span, y0, yp0, atol=atol, rtol=rtol, method=method, t_eval=t_eval)
end = time.time()
print(f"elapsed time: {end - start}")
t = sol.t
y = sol.y
tp = t[1:]
yp = np.diff(y) / np.diff(t)
yp = sol.yp
success = sol.success
status = sol.status
message = sol.message
Expand All @@ -96,15 +96,13 @@ def sol_true(t):
print(f"njev: {sol.njev}")
print(f"nlu: {sol.nlu}")

y_true, _ = sol_true(t)
_, yp_true = sol_true(tp)
y_true, yp_true = sol_true(t)
e = compute_error(y, y_true, rtol, atol)
ep = compute_error(yp, yp_true, rtol, atol)
# assert np.all(e < 5)
# print(f"e: {e}")
print(f"max(e): {max(e)}")
# print(f"max(ep): {max(ep)}")
# exit()

# visualization
fig, ax = plt.subplots(4, 1)
Expand All @@ -123,12 +121,12 @@ def sol_true(t):
ax[1].legend()
ax[1].grid()

ax[2].plot(tp, yp[4], "-ok", label="mup")
ax[2].plot(t, yp[4], "-ok", label="mup")
ax[2].legend()
ax[2].grid()

ax[3].plot(tp, yp[5], "-ok", label="lap")
ax[3].plot(tp, yp_true[5], "--xr", label="lap")
ax[3].plot(t, yp[5], "-ok", label="lap")
ax[3].plot(t, yp_true[5], "--xr", label="lap")
ax[3].legend()
ax[3].grid()

Expand Down
73 changes: 50 additions & 23 deletions examples/robertson.py
Original file line number Diff line number Diff line change
Expand Up @@ -10,41 +10,54 @@
References:
-----------
mathworks: https://de.mathworks.com/help/matlab/math/solve-differential-algebraic-equations-daes.html#bu75a7z-5 \\
Shampine2005: https://doi.org/10.1016/j.amc.2004.12.011
Shampine2005: https://doi.org/10.1016/j.amc.2004.12.011 \\
Sundials IDA (page 6): https://computing.llnl.gov/sites/default/files/ida_examples-5.7.0.pdf
"""

def f(t, y):
y1, y2, y3 = y

yp = np.zeros(3, dtype=y.dtype)
yp = np.zeros(3, dtype=float)
yp[0] = -0.04 * y1 + 1e4 * y2 * y3
yp[1] = 0.04 * y1 - 1e4 * y2 * y3 - 3e7 * y2**2
yp[2] = 3e7 * y2**2

return yp

def F(t, y, yp):
# return yp - f(t, y)
y1, y2, y3 = y
y1p, y2p, y3p = yp

F = np.zeros(3, dtype=y.dtype)
F = np.zeros(3, dtype=float)
F[0] = y1p - (-0.04 * y1 + 1e4 * y2 * y3)
F[1] = y2p - (0.04 * y1 - 1e4 * y2 * y3 - 3e7 * y2**2)
F[2] = y1 + y2 + y3 - 1

return F

def jac(t, y, yp):
y1, y2, y3 = y

Jyp = np.diag([1, 1, 0])
Jy = np.array([
[0.04, -1e4 * y3, -1e4 * y2],
[-0.04, 1e4 * y3 + 6e7 * y2, 1e4 * y2],
[1, 1, 1],
])
return Jy, Jyp


if __name__ == "__main__":
# time span
t0 = 0
t1 = 1e7
t1 = 4e10
t_span = (t0, t1)
t_eval = np.logspace(-6, 7, num=1000)
t_eval = np.logspace(-6, 7, num=200)
# t_eval = None

# method = "BDF"
method = "Radau"
method = "BDF"
# method = "Radau"

# initial conditions
y0 = np.array([1, 0, 0], dtype=float)
Expand All @@ -59,7 +72,7 @@ def F(t, y, yp):
print(f"fnorm: {fnorm}")

# solver options
atol = rtol = 1e-6
atol = rtol = 1e-8

####################
# reference solution
Expand All @@ -70,6 +83,7 @@ def F(t, y, yp):
print(f"elapsed time: {end - start}")
t_scipy = sol.t
y_scipy = sol.y
yp_scipy = np.array([f(ti, yi) for ti, yi in zip(t_scipy, y_scipy.T)]).T
success = sol.success
status = sol.status
message = sol.message
Expand All @@ -83,12 +97,14 @@ def F(t, y, yp):
##############
# dae solution
##############
stages = 3
start = time.time()
sol = solve_dae(F, t_span, y0, yp0, atol=atol, rtol=rtol, method=method, t_eval=t_eval)
sol = solve_dae(F, t_span, y0, yp0, atol=atol, rtol=rtol, method=method, t_eval=t_eval, jac=jac, stages=stages)
end = time.time()
print(f"elapsed time: {end - start}")
t = sol.t
y = sol.y
yp = sol.yp
success = sol.success
status = sol.status
message = sol.message
Expand All @@ -100,17 +116,28 @@ def F(t, y, yp):
print(f"nlu: {sol.nlu}")

# visualization
fig, ax = plt.subplots()

ax.set_xlabel("t")
ax.plot(t, y[0], "-ok", label="y1 DAE" + f" ({method})", mfc="none")
ax.plot(t, y[1] * 1e4, "-ob", label="y2 DAE" + f" ({method})", mfc="none")
ax.plot(t, y[2], "-og", label="y3 DAE" + f" ({method})", mfc="none")
ax.plot(t_scipy, y_scipy[0], "xr", label="y1 scipy Radau", markersize=7)
ax.plot(t_scipy, y_scipy[1] * 1e4, "xy", label="y2 scipy Radau", markersize=7)
ax.plot(t_scipy, y_scipy[2], "xm", label="y3 scipy Radau", markersize=7)
ax.set_xscale("log")
ax.legend()
ax.grid()

plt.show()
fig, ax = plt.subplots(2, 1)

ax[0].plot(t, y[0], "-ok", label="y1 DAE" + f" ({method})", mfc="none")
ax[0].plot(t, y[1] * 1e4, "-ob", label="y2 DAE" + f" ({method})", mfc="none")
ax[0].plot(t, y[2], "-og", label="y3 DAE" + f" ({method})", mfc="none")
ax[0].plot(t_scipy, y_scipy[0], "xr", label="y1 scipy" + f" ({method})", markersize=7)
ax[0].plot(t_scipy, y_scipy[1] * 1e4, "xy", label="y2 scipy" + f" ({method})", markersize=7)
ax[0].plot(t_scipy, y_scipy[2], "xm", label="y3 scipy" + f" ({method})", markersize=7)
ax[0].set_xlabel("t")
ax[0].set_xscale("log")
ax[0].legend()
ax[0].grid()

ax[1].plot(t, yp[0], "-ok", label="yp1 DAE" + f" ({method})", mfc="none")
ax[1].plot(t, yp[1], "-ob", label="yp2 DAE" + f" ({method})", mfc="none")
ax[1].plot(t, yp[2], "-og", label="yp3 DAE" + f" ({method})", mfc="none")
ax[1].plot(t_scipy, yp_scipy[0], "xr", label="yp1 scipy" + f" ({method})", markersize=7)
ax[1].plot(t_scipy, yp_scipy[1], "xy", label="yp2 scipy" + f" ({method})", markersize=7)
ax[1].plot(t_scipy, yp_scipy[2], "xm", label="yp3 scipy" + f" ({method})", markersize=7)
ax[1].set_xlabel("t")
ax[1].set_xscale("log")
ax[1].legend()
ax[1].grid()

plt.show()
65 changes: 63 additions & 2 deletions scipy_dae/integrate/_dae/base.py
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,6 @@
validate_first_step,
)
from scipy.optimize._numdiff import approx_derivative
from scipy.integrate._ivp.base import ConstantDenseOutput
from scipy.linalg import lu_factor, lu_solve
from scipy.sparse import csc_matrix, issparse, eye
from scipy.sparse.linalg import splu
Expand Down Expand Up @@ -419,7 +418,7 @@ def dense_output(self):

if self.n == 0 or self.t == self.t_old:
# Handle corner cases of empty solver and no integration.
return ConstantDenseOutput(self.t_old, self.t, self.y)
return ConstantDAEDenseOutput(self.t_old, self.t, self.y, self.yp)
else:
return self._dense_output_impl()

Expand All @@ -428,3 +427,65 @@ def _step_impl(self):

def _dense_output_impl(self):
raise NotImplementedError

class DAEDenseOutput:
"""Base class for local interpolant over step made by an ODE solver.
It interpolates between `t_min` and `t_max` (see Attributes below).
Evaluation outside this interval is not forbidden, but the accuracy is not
guaranteed.
Attributes
----------
t_min, t_max : float
Time range of the interpolation.
"""
def __init__(self, t_old, t):
self.t_old = t_old
self.t = t
self.t_min = min(t, t_old)
self.t_max = max(t, t_old)

def __call__(self, t):
"""Evaluate the interpolant.
Parameters
----------
t : float or array_like with shape (n_points,)
Points to evaluate the solution at.
Returns
-------
y : ndarray, shape (n,) or (n, n_points)
Computed values. Shape depends on whether `t` was a scalar or a
1-D array.
"""
t = np.asarray(t)
if t.ndim > 1:
raise ValueError("`t` must be a float or a 1-D array.")
return self._call_impl(t)

def _call_impl(self, t):
raise NotImplementedError


class ConstantDAEDenseOutput(DAEDenseOutput):
"""Constant value interpolator.
This class used for degenerate integration cases: equal integration limits
or a system with 0 equations.
"""
def __init__(self, t_old, t, value, derivative):
super().__init__(t_old, t)
self.value = value
self.derivative = derivative

def _call_impl(self, t):
if t.ndim == 0:
return self.value, self.derivative
else:
ret = np.empty((self.value.shape[0], t.shape[0]))
ret[:] = self.value[:, None]
ret_der = np.empty((self.derivative.shape[0], t.shape[0]))
ret_der[:] = self.derivative[:, None]
return ret, ret_der
91 changes: 75 additions & 16 deletions scipy_dae/integrate/_dae/bdf.py
Original file line number Diff line number Diff line change
@@ -1,7 +1,8 @@
import numpy as np
from math import factorial
from warnings import warn
from scipy.integrate._ivp.common import norm, EPS, warn_extraneous
from scipy.integrate._ivp.base import DenseOutput
from .base import DAEDenseOutput
from .dae import DaeSolver


Expand Down Expand Up @@ -220,6 +221,9 @@ class BDFDAE(DaeSolver):
.. [5] R. W. Klopfenstein, "Numerical differentiation formulas for stiff
systems of ordinary differential equations", RCA Review, 32,
pp. 447-462, September 1971.
.. [6] K. Radhakrishnan and A C. Hindmarsh, "Description and Use of LSODE,
the Livermore Solver for Ordinary Differential Equations", NASA
Reference Publication, December, 1993.
"""
def __init__(self, fun, t0, y0, yp0, t_bound, max_step=np.inf,
rtol=1e-3, atol=1e-6, jac=None, jac_sparsity=None,
Expand Down Expand Up @@ -418,26 +422,81 @@ def _dense_output_impl(self):
self.order, self.D[:self.order + 1].copy())


class BdfDenseOutput(DenseOutput):
class BdfDenseOutput(DAEDenseOutput):
def __init__(self, t_old, t, h, order, D):
super().__init__(t_old, t)
self.order = order
self.t_shift = self.t - h * np.arange(self.order)
self.denom = h * (1 + np.arange(self.order))
self.D = D
self.h = h

def _call_impl(self, t):
if t.ndim == 0:
x = (t - self.t_shift) / self.denom
p = np.cumprod(x)
else:
x = (t - self.t_shift[:, None]) / self.denom[:, None]
p = np.cumprod(x, axis=0)

y = np.dot(self.D[1:].T, p)
if y.ndim == 1:
y += self.D[0]
else:
y += self.D[0, :, None]

return y
vector_valued = t.ndim > 0
t = np.atleast_1d(t)

######################################################
# 0. default interpolation for y and yp is set to zero
######################################################
# x = (t - self.t_shift[:, None]) / self.denom[:, None]
# p = np.cumprod(x, axis=0)
# y = self.D[0, :, None] + np.dot(self.D[1:].T, p)

# # # TODO: Compute this derivative efficiently,
# # # see https://stackoverflow.com/questions/40916955/how-to-compute-gradient-of-cumprod-safely
# # dp = np.cumsum((p)[::-1], axis=0)[::-1] / x
# # yp = np.dot(self.D[1:].T, dp)
# yp = np.zeros_like(y)

# if vector_valued:
# y = np.squeeze(y)
# yp = np.squeeze(yp)

# return y, yp

############################################################
# 1. naive implementation of P(t) and P'(t) of p. 7 in [2]_.
############################################################
y2 = np.zeros((self.D.shape[1], len(t)), dtype=self.D.dtype)
y2 += self.D[0, :, None]
yp2 = np.zeros_like(y2)
for j in range(1, self.order + 1):
fac2 = np.ones_like(t)
dfac2 = np.zeros_like(t)
for m in range(j):
fac2 *= t - (self.t - m * self.h)

dprod2 = np.ones_like(t)
for i in range(j):
if i != m:
dprod2 *= t - (self.t - i * self.h)
dfac2 += dprod2

denom = factorial(j) * self.h**j
y2 += self.D[j, :, None] * fac2 / denom
yp2 += self.D[j, :, None] * dfac2 / denom

if vector_valued == 0:
y2 = np.squeeze(y2)
yp2 = np.squeeze(yp2)

return y2, yp2

#########################################################################
# 3. compute both values by Horner's rule,
# see see https://orionquest.github.io/Numacom/lectures/interpolation.pdf
#########################################################################
# y3 = np.zeros((self.D.shape[1], len(t)), dtype=self.D.dtype)
# y3 += self.D[n, :, None]
# yp3 = np.zeros_like(y3)
# for j in range(n, 0, -1):
# x = (t - (self.t - (j - 1) * self.h)) / (j * self.h)
# yp3 = y3 + x * yp3 * (j * self.h)
# yp3 /= (j * self.h)
# y3 = self.D[j - 1, :, None] + x * y3

# if vector_valued == 0:
# y3 = np.squeeze(y3)
# yp3 = np.squeeze(yp3)

# return y3, yp3
Loading

0 comments on commit fa0c9ee

Please sign in to comment.