From b402a28784c7f7bc287e3e2f692e8276152bde4d Mon Sep 17 00:00:00 2001 From: Jonas Breuling Date: Wed, 31 Jul 2024 11:02:07 +0200 Subject: [PATCH] Still no success with dense output for BDF and Robertson problem. --- examples/newton_polynomial.py | 3 +- examples/particle_on_circular_track.py | 4 +- examples/robertson.py | 4 +- scipy_dae/integrate/_dae/bdf.py | 102 ++++++++++++++++++++----- 4 files changed, 90 insertions(+), 23 deletions(-) diff --git a/examples/newton_polynomial.py b/examples/newton_polynomial.py index df70128..920c208 100644 --- a/examples/newton_polynomial.py +++ b/examples/newton_polynomial.py @@ -44,7 +44,8 @@ def newton_poly(coef, x_data, x): # sample points # x = np.array([-1, 0.25, 1]) # x = np.array([-2, -1, 0.25, 1, 2.5]) -num = 20 +num = 5 +# num = 20 x = np.linspace(-1, 1, num=num) # evaluate function values diff --git a/examples/particle_on_circular_track.py b/examples/particle_on_circular_track.py index efbe0cc..c39291c 100644 --- a/examples/particle_on_circular_track.py +++ b/examples/particle_on_circular_track.py @@ -75,8 +75,8 @@ def sol_true(t): t_eval = np.linspace(t0, t1, num=int(5e2)) # t_eval = None - # method = "BDF" - method = "Radau" + method = "BDF" + # method = "Radau" # initial conditions y0, yp0 = sol_true(t0) diff --git a/examples/robertson.py b/examples/robertson.py index 3686b3c..4faabaa 100644 --- a/examples/robertson.py +++ b/examples/robertson.py @@ -59,8 +59,8 @@ def fun_composite(t, z): # t_eval = np.logspace(-6, 7, num=1000) t_eval = np.logspace(-6, 7, num=500) - # method = "BDF" - method = "Radau" + method = "BDF" + # method = "Radau" # initial conditions y0 = np.array([1, 0, 0], dtype=float) diff --git a/scipy_dae/integrate/_dae/bdf.py b/scipy_dae/integrate/_dae/bdf.py index eaf4b1a..7d0179b 100644 --- a/scipy_dae/integrate/_dae/bdf.py +++ b/scipy_dae/integrate/_dae/bdf.py @@ -365,6 +365,7 @@ def _step_impl(self): else: step_accepted = True + print(f"step accepted at t: {t} with order: {order}") self.n_equal_steps += 1 self.t = t_new @@ -418,6 +419,7 @@ def _step_impl(self): # delta_order = np.argmin(error_norms) - 1 order += delta_order + # order = max(2, order) # TODO: Remove this later self.order = order factor = min(MAX_FACTOR, safety * np.max(factors)) @@ -468,24 +470,86 @@ def _call_impl(self, t): yp2 = np.zeros_like(y) for j in range(1, self.order + 1): # interpolate y - factor_y2 = np.ones_like(t) - factor_yp2 = np.zeros_like(t) + fac2 = np.ones_like(t) + dfac2 = np.zeros_like(t) for m in range(j): - factor_y2 *= t - (self.t - m * self.h) - product2 = np.ones_like(t) + fac2 *= t - (self.t - m * self.h) + prod2 = np.ones_like(t) for i in range(j): if i != m: - product2 *= t - (self.t - m * self.h) - factor_yp2 += product2 + # prod2 *= t - (self.t - m * self.h) + prod2 *= t - (self.t - i * self.h) + dfac2 += prod2 if y.ndim == 1: - y2 += self.D[j] * factor_y2 / (factorial(j) * self.h**j) - yp2 += self.D[j] * factor_yp2 / (factorial(j) * self.h**j) + y2 += self.D[j] * fac2 / (factorial(j) * self.h**j) + yp2 += self.D[j] * dfac2 / (factorial(j) * self.h**j) else: - y2 += self.D[j, :, None] * factor_y2 / (factorial(j) * self.h**j) - yp2 += self.D[j, :, None] * factor_yp2 / (factorial(j) * self.h**j) + y2 += self.D[j, :, None] * fac2 / (factorial(j) * self.h**j) + yp2 += self.D[j, :, None] * dfac2 / (factorial(j) * self.h**j) assert np.allclose(y, y2) + # return y2, yp2 + + # Shampine p. 7 + y3 = np.zeros_like(y) + if y.ndim == 1: + y3 += self.D[0] + else: + y3 += self.D[0, :, None] + yp3 = np.zeros_like(y) + for j in range(1, self.order + 1): + fac3 = np.ones_like(t) + dfac3 = np.zeros_like(t) + for m in range(j): + fac3 *= t - (self.t - m * self.h) + + dprod3 = np.ones_like(t) + for i in range(j): + if i != m: + # dprod3 *= t - (self.t - m * self.h) + dprod3 *= t - (self.t - i * self.h) + dfac3 += dprod3 + + fac3 /= factorial(j) * self.h**j + dfac3 /= factorial(j) * self.h**j + + if y.ndim == 1: + y3 += self.D[j] * fac3 + yp3 += self.D[j] * dfac3 + else: + y3 += self.D[j, :, None] * fac3[None, :] + yp3 += self.D[j, :, None] * dfac3[None, :] + + assert np.allclose(y, y3) + + # return y3, yp3 + + # Horners method, see https://orionquest.github.io/Numacom/lectures/interpolation.pdf + n = self.order + y4 = np.zeros_like(y) + if y.ndim == 1: + y4 += self.D[n] + else: + y4 += self.D[n, :, None] + yp4 = np.zeros_like(y) + + # for j in range(self.order, -1, -1): + # for j in range(n, -1, -1): + for j in range(n - 1, -1, -1): + x = (t - (self.t - j * self.h)) / ((j + 1) * self.h) + # x = (t - (self.t - j * self.h)) / self.h + # x = (t - (self.t - j * self.h))# / self.h + # x = (t - self.t) #/ self.h + yp4 = y4 + x * yp4 + yp4 /= ((j + 1) * self.h) + if y.ndim == 1: + y4 = self.D[j] + x * y4 + else: + y4 = self.D[j, :, None] + x * y4 + # yp4 /= factorial(n - 1) * self.h#**(n - 1) + # yp4 /= self.h + # return y4, yp4 # compute y recursively which enables the computation of yp as well using Horner's rule, see # https://pages.cs.wisc.edu/~amos/412/lecture-notes/lecture08.pdf @@ -500,14 +564,16 @@ def _call_impl(self, t): yp3 = y3 * factor + p * yp3 # TODO: Understand factor * y3 here y3 = self.D[j - 1, :, None] + p * y3 - # print(f"|y - y2|: {np.linalg.norm(y - y2)}") - # # print(f"|yp - yp2|: {np.linalg.norm(yp - yp2)}") + return y3, yp3 #/ self.h - # print(f"|y - y3|: {np.linalg.norm(y - y3)}") - print(f"|yp2 - yp3|: {np.linalg.norm(yp2 - yp3)}") - assert np.allclose(y, y3) - # assert np.allclose(yp2, yp3) - y = y2 - yp = yp2 + # # print(f"|y - y2|: {np.linalg.norm(y - y2)}") + # # # print(f"|yp - yp2|: {np.linalg.norm(yp - yp2)}") + + # # print(f"|y - y3|: {np.linalg.norm(y - y3)}") + # print(f"|yp2 - yp3|: {np.linalg.norm(yp2 - yp3)}") + # assert np.allclose(y, y3) + # # assert np.allclose(yp2, yp3) + # y = y2 + # yp = yp2 return y, yp