From 8518e9c84b5533df1339b5083d5f20aff445726b Mon Sep 17 00:00:00 2001 From: Jonas Breuling Date: Tue, 30 Jul 2024 15:29:41 +0200 Subject: [PATCH] Derivative of interpolation polynomial is also correct. --- scipy_dae/integrate/_dae/radau.py | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/scipy_dae/integrate/_dae/radau.py b/scipy_dae/integrate/_dae/radau.py index dcebbee..c450484 100644 --- a/scipy_dae/integrate/_dae/radau.py +++ b/scipy_dae/integrate/_dae/radau.py @@ -699,8 +699,14 @@ def __init__(self, t_old, t, y_old, Q, yp_old, Qp): 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) + + # factors for nterpolation polynomial and its derivative + # p = np.tile(x, (self.order + 1, 1)) + # p = np.cumprod(p, axis=0) + c = np.arange(1, self.order + 2)[:, None] + p = x**c + dp = (c / self.h) * (x**(c - 1)) + # Here we don't multiply by h, not a mistake. y = np.dot(self.Q, p) yp = np.dot(self.Qp, p) @@ -710,6 +716,14 @@ def _call_impl(self, t): y = np.squeeze(y) yp = np.squeeze(yp) + # compute derivative of interpolation polynomial + yp = np.dot(self.Q, dp) + + # # compute derivative of interpolation polynomial + # yp = np.zeros_like(y) + # for i in range(self.order + 1): + # yp += np.outer(self.Q[:, i] * (i + 1), x**i) / self.h + return y, yp