From fa0c9ee6a010402115a19cb5efa7ec9c3836f266 Mon Sep 17 00:00:00 2001 From: Jonas Breuling Date: Fri, 2 Aug 2024 14:40:57 +0200 Subject: [PATCH] Dense output for yp. --- examples/particle_on_circular_track.py | 18 ++- examples/robertson.py | 73 +++++++++---- scipy_dae/integrate/_dae/base.py | 65 ++++++++++- scipy_dae/integrate/_dae/bdf.py | 91 +++++++++++++--- scipy_dae/integrate/_dae/common.py | 121 +++++++++++++++++++++ scipy_dae/integrate/_dae/dae.py | 41 +++---- scipy_dae/integrate/_dae/radau.py | 46 ++++++-- scipy_dae/integrate/_dae/tests/test_dae.py | 6 +- 8 files changed, 373 insertions(+), 88 deletions(-) diff --git a/examples/particle_on_circular_track.py b/examples/particle_on_circular_track.py index 8231a50..112b07c 100644 --- a/examples/particle_on_circular_track.py +++ b/examples/particle_on_circular_track.py @@ -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" @@ -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 @@ -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) @@ -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() diff --git a/examples/robertson.py b/examples/robertson.py index d2c562a..0f54d90 100644 --- a/examples/robertson.py +++ b/examples/robertson.py @@ -10,13 +10,14 @@ 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 @@ -24,27 +25,39 @@ def f(t, y): 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) @@ -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 @@ -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 @@ -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 @@ -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() \ No newline at end of file diff --git a/scipy_dae/integrate/_dae/base.py b/scipy_dae/integrate/_dae/base.py index fc99025..284b426 100644 --- a/scipy_dae/integrate/_dae/base.py +++ b/scipy_dae/integrate/_dae/base.py @@ -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 @@ -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() @@ -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 diff --git a/scipy_dae/integrate/_dae/bdf.py b/scipy_dae/integrate/_dae/bdf.py index 20a83b4..de3ae54 100644 --- a/scipy_dae/integrate/_dae/bdf.py +++ b/scipy_dae/integrate/_dae/bdf.py @@ -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 @@ -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, @@ -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 diff --git a/scipy_dae/integrate/_dae/common.py b/scipy_dae/integrate/_dae/common.py index 7026ccd..78d637a 100644 --- a/scipy_dae/integrate/_dae/common.py +++ b/scipy_dae/integrate/_dae/common.py @@ -1,3 +1,4 @@ +from itertools import groupby import numpy as np # TODO: use sparse QR when available in scipy from scipy.linalg import qr, solve_triangular @@ -5,6 +6,126 @@ from scipy.optimize._numdiff import approx_derivative +class DaeSolution: + """Continuous DAE solution. + It is organized as a collection of `DenseOutput` objects which represent + local interpolants. It provides an algorithm to select a right interpolant + for each given point. + The interpolants cover the range between `t_min` and `t_max` (see + Attributes below). Evaluation outside this interval is not forbidden, but + the accuracy is not guaranteed. + When evaluating at a breakpoint (one of the values in `ts`) a segment with + the lower index is selected. + Parameters + ---------- + ts : array_like, shape (n_segments + 1,) + Time instants between which local interpolants are defined. Must + be strictly increasing or decreasing (zero segment with two points is + also allowed). + interpolants : list of DenseOutput with n_segments elements + Local interpolants. An i-th interpolant is assumed to be defined + between ``ts[i]`` and ``ts[i + 1]``. + alt_segment : boolean + Requests the alternative interpolant segment selection scheme. At each + solver integration point, two interpolant segments are available. The + default (False) and alternative (True) behaviours select the segment + for which the requested time corresponded to ``t`` and ``t_old``, + respectively. This functionality is only relevant for testing the + interpolants' accuracy: different integrators use different + construction strategies. + Attributes + ---------- + t_min, t_max : float + Time range of the interpolation. + """ + def __init__(self, ts, interpolants, alt_segment=False): + ts = np.asarray(ts) + d = np.diff(ts) + # The first case covers integration on zero segment. + if not ((ts.size == 2 and ts[0] == ts[-1]) + or np.all(d > 0) or np.all(d < 0)): + raise ValueError("`ts` must be strictly increasing or decreasing.") + + self.n_segments = len(interpolants) + if ts.shape != (self.n_segments + 1,): + raise ValueError("Numbers of time stamps and interpolants " + "don't match.") + + self.ts = ts + self.interpolants = interpolants + if ts[-1] >= ts[0]: + self.t_min = ts[0] + self.t_max = ts[-1] + self.ascending = True + self.side = "right" if alt_segment else "left" + self.ts_sorted = ts + else: + self.t_min = ts[-1] + self.t_max = ts[0] + self.ascending = False + self.side = "left" if alt_segment else "right" + self.ts_sorted = ts[::-1] + + def _call_single(self, t): + # Here we preserve a certain symmetry that when t is in self.ts, + # if alt_segment=False, then we prioritize a segment with a lower + # index. + ind = np.searchsorted(self.ts_sorted, t, side=self.side) + + segment = min(max(ind - 1, 0), self.n_segments - 1) + if not self.ascending: + segment = self.n_segments - 1 - segment + + return self.interpolants[segment](t) + + def __call__(self, t): + """Evaluate the solution. + Parameters + ---------- + t : float or array_like with shape (n_points,) + Points to evaluate at. + Returns + ------- + y : ndarray, shape (n_states,) or (n_states, n_points) + Computed values. Shape depends on whether `t` is a scalar or a + 1-D array. + """ + t = np.asarray(t) + + if t.ndim == 0: + return self._call_single(t) + + order = np.argsort(t) + reverse = np.empty_like(order) + reverse[order] = np.arange(order.shape[0]) + t_sorted = t[order] + + # See comment in self._call_single. + segments = np.searchsorted(self.ts_sorted, t_sorted, side=self.side) + segments -= 1 + segments[segments < 0] = 0 + segments[segments > self.n_segments - 1] = self.n_segments - 1 + if not self.ascending: + segments = self.n_segments - 1 - segments + + ys = [] + yps = [] + group_start = 0 + for segment, group in groupby(segments): + group_end = group_start + len(list(group)) + y, yp = self.interpolants[segment](t_sorted[group_start:group_end]) + ys.append(y) + yps.append(yp) + group_start = group_end + + ys = np.hstack(ys) + ys = ys[:, reverse] + yps = np.hstack(yps) + yps = yps[:, reverse] + + return ys, yps + + # TODO: Compare this with # - ddassl.f by Petzold # - epsode.f by Bryne and Hindmarsh diff --git a/scipy_dae/integrate/_dae/dae.py b/scipy_dae/integrate/_dae/dae.py index cefb732..c46f233 100644 --- a/scipy_dae/integrate/_dae/dae.py +++ b/scipy_dae/integrate/_dae/dae.py @@ -5,11 +5,12 @@ import inspect import numpy as np -from scipy.integrate._ivp.ivp import OdeResult, prepare_events, handle_events, find_active_events -from scipy.integrate._ivp.common import OdeSolution +from scipy.integrate._ivp.ivp import prepare_events, handle_events, find_active_events from .base import DaeSolver from .bdf import BDFDAE from .radau import RadauDAE +from scipy.optimize import OptimizeResult +from .common import DaeSolution METHODS = { @@ -22,11 +23,12 @@ 1: "A termination event occurred."} +class DaeResult(OptimizeResult): + pass + + # TODO: -# - expect consistent initial conditions and add a helper function that -# computes them as done by matlab? # - add events depending on y'(t)? -# - dense output for y' def solve_dae(fun, t_span, y0, y_dot0, method="Radau", t_eval=None, dense_output=False, events=None, vectorized=False, args=None, **options): @@ -203,8 +205,8 @@ def solve_dae(fun, t_span, y0, y_dot0, method="Radau", t_eval=None, Time points. y : ndarray, shape (n, n_points) Values of the solution at `t`. - sol : `OdeSolution` or None - Found solution as `OdeSolution` instance; None if `dense_output` was + sol : `DaeSolution` or None + Found solution as `DaeSolution` instance; None if `dense_output` was set to False. t_events : list of ndarray or None Contains for each event type a list of arrays at which an event of @@ -413,7 +415,9 @@ def fun(t, x, x_dot, fun=fun): if sol is None: sol = solver.dense_output() ts.append(t_eval_step) - ys.append(sol(t_eval_step)) + y, yp = sol(t_eval_step) + ys.append(y) + yps.append(yp) t_eval_i = t_eval_i_new if t_eval is not None and dense_output: @@ -433,29 +437,16 @@ def fun(t, x, x_dot, fun=fun): elif ts: ts = np.hstack(ts) ys = np.hstack(ys) - # yps = np.hstack(yps) + yps = np.hstack(yps) - # estimate yps via finite differences - t_final1 = 2 * ts[-1] - ts[-2] - y_final1 = 2 * ys[:, -1] - ys[:, -2] - ts1 = np.array([*ts, t_final1]) - ys1 = np.concatenate((ys, y_final1[:, None]), axis=1) - yps = np.diff(ys1) / np.diff(ts1) - if dense_output: if t_eval is None: - sol = OdeSolution( - ts, interpolants, #alt_segment=False - # ts, interpolants, alt_segment=True if method in [BDFDAE] else False - ) + sol = DaeSolution(ts, interpolants) else: - sol = OdeSolution( - ti, interpolants, #alt_segment=False - # ti, interpolants, alt_segment=True if method in [BDFDAE] else False - ) + sol = DaeSolution(ti, interpolants) else: sol = None - return OdeResult(t=ts, y=ys, yp=yps, sol=sol, t_events=t_events, y_events=y_events, + return DaeResult(t=ts, y=ys, yp=yps, sol=sol, t_events=t_events, y_events=y_events, nfev=solver.nfev, njev=solver.njev, nlu=solver.nlu, status=status, message=message, success=status>=0) diff --git a/scipy_dae/integrate/_dae/radau.py b/scipy_dae/integrate/_dae/radau.py index cf44d60..6d5f8ea 100644 --- a/scipy_dae/integrate/_dae/radau.py +++ b/scipy_dae/integrate/_dae/radau.py @@ -2,7 +2,7 @@ from numpy.polynomial import Polynomial as Poly from scipy.linalg import eig, cdf2rdf 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 @@ -499,7 +499,7 @@ def _step_impl(self): if self.sol is None: Z0 = np.zeros((s, y.shape[0])) else: - Z0 = self.sol(t + h * C).T - y + Z0 = self.sol(t + h * C)[0].T - y scale = atol + np.abs(y) * rtol converged = False @@ -624,29 +624,57 @@ def _step_impl(self): def _compute_dense_output(self): Q = np.dot(self.Z.T, self.P) - return RadauDenseOutput(self.t_old, self.t, self.y_old, Q) + h = self.t - self.t_old + Yp = (self.A_inv / h) @ self.Z + Zp = Yp - self.yp_old + Qp = np.dot(Zp.T, self.P) + return RadauDenseOutput(self.t_old, self.t, self.y_old, Q, self.yp_old, Qp) def _dense_output_impl(self): return self.sol -class RadauDenseOutput(DenseOutput): - def __init__(self, t_old, t, y_old, Q): +class RadauDenseOutput(DAEDenseOutput): + def __init__(self, t_old, t, y_old, Q, yp_old, Qp): super().__init__(t_old, t) self.h = t - t_old self.Q = Q + self.Qp = Qp self.order = Q.shape[1] - 1 self.y_old = y_old + self.yp_old = yp_old def _call_impl(self, t): x = (t - self.t_old) / self.h x = np.atleast_1d(x) - p = np.tile(x, (self.order + 1, 1)) - p = np.cumprod(p, axis=0) - # Here we don't multiply by h, not a mistake. + + # factors for interpolation polynomial and its derivative + c = np.arange(1, self.order + 2)[:, None] + p = x**c + dp = (c / self.h) * (x**(c - 1)) + + # # 1. compute derivative of interpolation polynomial for y y = np.dot(self.Q, p) y += self.y_old[:, None] + yp = np.dot(self.Q, dp) + + # # 2. compute collocation polynomial for y and yp + # y = np.dot(self.Q, p) + # yp = np.dot(self.Qp, p) + # y += self.y_old[:, None] + # yp += self.yp_old[:, None] + + # # 3. compute both values by Horner's rule + # y = np.zeros_like(y) + # yp = np.zeros_like(y) + # for i in range(self.order, -1, -1): + # y = self.Q[:, i][:, None] + y * x[None, :] + # yp = y + yp * x[None, :] + # y = self.y_old[:, None] + y * x[None, :] + # yp /= self.h + if t.ndim == 0: y = np.squeeze(y) + yp = np.squeeze(yp) - return y + return y, yp diff --git a/scipy_dae/integrate/_dae/tests/test_dae.py b/scipy_dae/integrate/_dae/tests/test_dae.py index 88ebc6b..202827c 100644 --- a/scipy_dae/integrate/_dae/tests/test_dae.py +++ b/scipy_dae/integrate/_dae/tests/test_dae.py @@ -174,19 +174,19 @@ def test_integration_rational(vectorized, method, t_span, jac): tc = np.linspace(*t_span) yc_true = sol_rational(tc) - yc = res.sol(tc) + yc = res.sol(tc)[0] e = compute_error(yc, yc_true, rtol, atol) assert_(np.all(e < 6)) tc = (t_span[0] + t_span[-1]) / 2 yc_true = sol_rational(tc) - yc = res.sol(tc) + yc = res.sol(tc)[0] e = compute_error(yc, yc_true, rtol, atol) assert_(np.all(e < 5)) - assert_allclose(res.sol(res.t), res.y, rtol=1e-15, atol=1e-15) + assert_allclose(res.sol(res.t)[0], res.y, rtol=1e-15, atol=1e-15) parameters_stiff = ["BDF", "Radau"]