Skip to content

Commit

Permalink
Cleanup Radau again.
Browse files Browse the repository at this point in the history
  • Loading branch information
@JonasBreuling
JonasBreuling committed Aug 15, 2024
1 parent 80ee92e commit 12b3c10
Show file tree
Hide file tree
Showing 2 changed files with 48 additions and 117 deletions.
0 examples/daes/one_element_timoshenko_rod.py → examples/odes/one_element_timoshenko_rod.py
View file
File renamed without changes.
165 changes: 48 additions & 117 deletions scipy_dae/integrate/_dae/radau.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -6,15 +6,11 @@
from .dae import DaeSolver


NEWTON_MAXITER_EMBEDDED = 3 # Maximum number of Newton iterations for embedded method.
DAMPING_RATIO_ERROR_ESTIMATE = 0.8 # Hairer (8.19) is obtained by the choice 1.0.
# de Swart proposes 0.067 for s=3.
MIN_FACTOR = 0.2 # Minimum allowed decrease in a step size.
MAX_FACTOR = 10 # Maximum allowed increase in a step size.
KAPPA = 1 # Factor of the smooth limiter
# VELOCITIY_FORMULATION = False
VELOCITIY_FORMULATION = True
SCALED_VELOCITIES = True


def butcher_tableau(s):
Expand Down Expand Up @@ -46,10 +42,7 @@ def radau_constants(s):
A_inv = np.linalg.inv(A)

# eigenvalues and corresponding eigenvectors of inverse coefficient matrix
if VELOCITIY_FORMULATION:
lambdas, V = eig(A)
else:
lambdas, V = eig(A_inv)
lambdas, V = eig(A_inv)

# sort eigenvalues and permut eigenvectors accordingly
idx = np.argsort(lambdas)[::-1]
Expand All @@ -67,16 +60,10 @@ def radau_constants(s):
TI = np.linalg.inv(T)

# check if everything worked
if VELOCITIY_FORMULATION:
assert np.allclose(V @ np.diag(lambdas) @ np.linalg.inv(V), A)
assert np.allclose(np.linalg.inv(V) @ A @ V, np.diag(lambdas))
assert np.allclose(T @ Gammas @ TI, A)
assert np.allclose(TI @ A @ T, Gammas)
else:
assert np.allclose(V @ np.diag(lambdas) @ np.linalg.inv(V), A_inv)
assert np.allclose(np.linalg.inv(V) @ A_inv @ V, np.diag(lambdas))
assert np.allclose(T @ Gammas @ TI, A_inv)
assert np.allclose(TI @ A_inv @ T, Gammas)
assert np.allclose(V @ np.diag(lambdas) @ np.linalg.inv(V), A_inv)
assert np.allclose(np.linalg.inv(V) @ A_inv @ V, np.diag(lambdas))
assert np.allclose(T @ Gammas @ TI, A_inv)
assert np.allclose(TI @ A_inv @ T, Gammas)

# extract real and complex eigenvalues
real_idx = lambdas.imag == 0
Expand All @@ -91,11 +78,7 @@ def radau_constants(s):
vander = np.vander(c_hat, increasing=True).T

rhs = 1 / np.arange(1, s + 1)
if VELOCITIY_FORMULATION:
# TODO:
b_hats2 = 1 / gammas[0] # real eigenvalue of A, i.e., 1 / gamma[0]
else:
b_hats2 = 1 / gammas[0] # real eigenvalue of A, i.e., 1 / gamma[0]
b_hats2 = 1 / gammas[0] # real eigenvalue of A, i.e., 1 / gamma[0]
b_hat1 = DAMPING_RATIO_ERROR_ESTIMATE * b_hats2
rhs[0] -= b_hat1
rhs -= b_hats2
Expand Down Expand Up @@ -175,21 +158,10 @@ def solve_collocation_system(fun, t, y, h, Z0, scale, tol,
ncs = s // 2
tau = t + h * C

if VELOCITIY_FORMULATION:
if SCALED_VELOCITIES:
Yp = Z0
Z = Yp * h
W = TI.dot(Z)
Y = y + A @ Z
else:
Yp = Z0
W = TI.dot(Z0)
Y = y + h * A @ Yp
else:
Z = Z0
W = TI.dot(Z0)
Yp = (A_inv / h) @ Z
Y = y + Z
Z = Z0
W = TI.dot(Z0)
Yp = (A_inv / h) @ Z
Y = y + Z

F = np.empty((s, n))

Expand Down Expand Up @@ -244,18 +216,9 @@ def solve_collocation_system(fun, t, y, h, Z0, scale, tol,
# break

W += dW
if VELOCITIY_FORMULATION:
if SCALED_VELOCITIES:
Z = T.dot(W)
Yp = Z / h
Y = y + A @ Z
else:
Yp = T.dot(W)
Y = y + h * A @ Yp
else:
Z = T.dot(W)
Yp = (A_inv / h) @ Z
Y = y + Z
Z = T.dot(W)
Yp = (A_inv / h) @ Z
Y = y + Z

if (dW_norm == 0 or rate is not None and rate / (1 - rate) * dW_norm < tol):
converged = True
Expand All @@ -278,14 +241,6 @@ def solve_collocation_system(fun, t, y, h, Z0, scale, tol,
dW_norm_old = dW_norm
rate_old = rate

if VELOCITIY_FORMULATION:
if SCALED_VELOCITIES:
# Y = y + A @ Z
# Z = A_inv @ Z
Z = Y - y
else:
Z = Y - y

return converged, k + 1, Y, Yp, Z, rate


Expand Down Expand Up @@ -584,7 +539,7 @@ def __init__(self, fun, t0, y0, yp0, t_bound, stages=3,
continuous_error_weight=0.0, jac=None,
jac_sparsity=None, vectorized=False,
first_step=None, newton_max_iter=None,
jac_recompute_rate=1e-3,
jac_recompute_rate=1e-3, newton_iter_embedded=1,
controller_deadband=(1.0, 1.2),
**extraneous):
warn_extraneous(extraneous)
Expand Down Expand Up @@ -641,6 +596,9 @@ def __init__(self, fun, t0, y0, yp0, t_bound, stages=3,
assert 0 < controller_deadband[0] <= controller_deadband[1]
self.controller_deadband = controller_deadband

assert 0 <= newton_iter_embedded
self.newton_iter_embedded = newton_iter_embedded

self.current_jac = True
self.LU_real = None
self.LU_complex = None
Expand Down Expand Up @@ -713,41 +671,17 @@ def _step_impl(self):
h = t_new - t
h_abs = np.abs(h)

if VELOCITIY_FORMULATION:
if self.sol is None:
Z0 = np.zeros((s, y.shape[0]))
else:
# actually, this is Yp0
Z0 = self.sol(t + h * C)[1].T
# Z0 = np.zeros((s, y.shape[0]))
# Z0 = self.sol(t + h * C)[0].T - y
# Z0 = (A_inv / h) @ Z0
if SCALED_VELOCITIES:
# scale = atol + np.abs(h * yp) * rtol
scale = atol + np.abs(y) * rtol
else:
scale = atol + np.abs(yp) * rtol
if self.sol is None:
Z0 = np.zeros((s, y.shape[0]))
else:
if self.sol is None:
Z0 = np.zeros((s, y.shape[0]))
else:
Z0 = self.sol(t + h * C)[0].T - y
Z0 = np.zeros((s, y.shape[0]))
scale = atol + np.abs(y) * rtol
Z0 = self.sol(t + h * C)[0].T - y
scale = atol + np.abs(y) * rtol

converged = False
while not converged:
if LU_real is None or LU_complex is None:
if VELOCITIY_FORMULATION:
if SCALED_VELOCITIES:
LU_real = self.lu(MU_REAL * Jy + (1 / h) * Jyp)
LU_complex = [self.lu(MU * Jy + (1 / h) * Jyp) for MU in MU_COMPLEX]
else:
LU_real = self.lu(MU_REAL * h * Jy + Jyp)
LU_complex = [self.lu(MU * h * Jy + Jyp) for MU in MU_COMPLEX]
else:
LU_real = self.lu(MU_REAL / h * Jyp + Jy)
LU_complex = [self.lu(MU / h * Jyp + Jy) for MU in MU_COMPLEX]
LU_real = self.lu(MU_REAL / h * Jyp + Jy)
LU_complex = [self.lu(MU / h * Jyp + Jy) for MU in MU_COMPLEX]

converged, n_iter, Y, Yp, Z, rate = solve_collocation_system(
self.fun, t, y, h, Z0, scale, newton_tol,
Expand Down Expand Up @@ -775,34 +709,34 @@ def _step_impl(self):

scale = atol + np.maximum(np.abs(y), np.abs(y_new)) * rtol

############################################
# error of collocation polynomial of order s
############################################
error_collocation = y - P2[0] @ Y
error_embedded = error_collocation

# # explicit embedded method with R(z) = +oo for z -> oo
# error_embedded = h * (yp / MU_REAL + v2 @ Yp)
# embedded error measures
if self.newton_iter_embedded == 0:
# explicit embedded method with R(z) = +oo for z -> oo
error_embedded = h * (yp / MU_REAL + v2 @ Yp)
elif self.newton_iter_embedded == 1:

# # compute implicit embedded method with a single iteration
# # # TODO: Store MU_REAL * v during construction.
# # error_embedded = h * MU_REAL * (v @ Yp - b0 * yp - yp_new / MU_REAL)
# yp_hat_new = MU_REAL * (v @ Yp - b0 * yp)
# F = self.fun(t_new, y_new, yp_hat_new)
# error_embedded = self.solve_lu(LU_real, -F)

# # compute implicit embedded method with NEWTON_MAXITER_EMBEDDED iterations
# y_hat_new = y_new.copy() # initial guess
# for _ in range(NEWTON_MAXITER_EMBEDDED):
# yp_hat_new = MU_REAL * (
# (y_hat_new - y) / h
# - b0 * yp
# - self.b_hat @ Yp
# )
# F = self.fun(t_new, y_hat_new, yp_hat_new)
# y_hat_new -= self.solve_lu(LU_real, F)

# error_embedded = y_hat_new - y_new
# compute implicit embedded method with a single Newton iteration;
# R(z) = b_hat1 / b_hats2 = DAMPING_RATIO_ERROR_ESTIMATE for z -> oo
# TODO: Store MU_REAL * v during construction.
yp_hat_new = MU_REAL * (v @ Yp - b0 * yp)
F = self.fun(t_new, y_new, yp_hat_new)
error_embedded = self.solve_lu(LU_real, -F)
else:
# compute implicit embedded method with `newton_iter_embedded`` iterations
y_hat_new = y_new.copy() # initial guess
for _ in range(self.newton_iter_embedded):
yp_hat_new = MU_REAL * (
(y_hat_new - y) / h
- b0 * yp
- self.b_hat @ Yp
)
F = self.fun(t_new, y_hat_new, yp_hat_new)
y_hat_new -= self.solve_lu(LU_real, F)

error_embedded = y_hat_new - y_new

# mix embedded error with collocation error as proposed in Guglielmi2001/ Guglielmi2003
error = (
Expand Down Expand Up @@ -836,10 +770,7 @@ def _step_impl(self):
# factor = min(MAX_FACTOR, safety * factor)
factor = safety * factor

# TODO: We should also use the ideas of Söderlind
# (https://doi.org/10.1016/j.cam.2005.03.008)
# that proposes a smooth change of the step-size.
# if not recompute_jac and factor < 1.2:
# do not alter step-size in deadband
if (
not recompute_jac
and self.controller_deadband[0] <= factor <= self.controller_deadband[1]
Expand Down

0 comments on commit 12b3c10

Please sign in to comment.