diff --git a/README.md b/README.md index a7f6ae3..b900fca 100644 --- a/README.md +++ b/README.md @@ -64,8 +64,7 @@ EverhartIILobatto21ODESolver ____ ## Example: -```python -import math +```pythonimport math import numpy as np import matplotlib.pyplot as plt @@ -159,6 +158,7 @@ if __name__ == '__main__': dt0 = np.longdouble(0.01) tolerance = 1e-4 + is_adaptive_step = True # initial conditions: x0 = np.longdouble(0) @@ -169,7 +169,7 @@ if __name__ == '__main__': vz0 = np.longdouble(0) u0 = np.array([x0, vx0, y0, vy0, z0, vz0], dtype='longdouble') - solver = Fehlberg45Solver(f, u0, t0, tmax, dt0, is_adaptive_step=True, tolerance=tolerance) + solver = Fehlberg45Solver(f, u0, t0, tmax, dt0, is_adaptive_step=is_adaptive_step, tolerance=tolerance) solution, time_points = solver.solve(print_benchmark=True, benchmark_name=solver.name) x_solution = solution[:, 0] y_solution = solution[:, 2] @@ -183,7 +183,7 @@ if __name__ == '__main__': u0 = np.array([x0, y0, z0], dtype='longdouble') du_dt0 = np.array([vx0, vy0, vz0], dtype='longdouble') - solver1 = EverhartIIRadau7ODESolver(f2, u0, du_dt0, t0, tmax, dt0, is_adaptive_step=True, + solver1 = EverhartIIRadau7ODESolver(f2, u0, du_dt0, t0, tmax, dt0, is_adaptive_step=is_adaptive_step, tolerance=tolerance) solution1, d_solution1, time_points1 = solver1.solve(print_benchmark=True, benchmark_name=solver1.name) x_solution1 = solution1[:, 0] @@ -193,7 +193,12 @@ if __name__ == '__main__': ax.plot(time_points1, x_solution1, label=solver1.name) points_number = int((tmax - t0) / dt0) - time_exact = np.linspace(t0, t0 + dt0 * points_number, (points_number + 1)) + + if is_adaptive_step: + time_exact = np.linspace(t0, t0 + dt0 * points_number, (points_number + 1)) + else: + time_exact = np.linspace(t0, t0 + dt0 * points_number, points_number + 2) + x_exact, y_exact = exact_f(time_exact) plt.plot(time_exact, x_exact, label='Exact analytical solution') diff --git a/cleanode/ode_solvers.py b/cleanode/ode_solvers.py index f971433..2d7fda9 100644 --- a/cleanode/ode_solvers.py +++ b/cleanode/ode_solvers.py @@ -99,6 +99,7 @@ def solve(self, print_benchmark=False, benchmark_name='') -> Tuple[np.ndarray, n if self.is_adaptive_step: u_next = None real_tolerance = self.tolerance * 2 + dt_current = self.dt while real_tolerance > self.tolerance: dt_current = self.dt # noinspection PyTypeChecker @@ -874,7 +875,7 @@ def __init__(self, order, u0 = np.asarray(u0, dtype='longdouble') self.ode_system_size = u0.size - self.alfa = np.zeros([len(self.h), len(u0)], dtype='longdouble') + self.alfa = np.zeros([len(self.h) - 1, len(u0)], dtype='longdouble') self.dt = dt0 @@ -910,7 +911,6 @@ def solve(self, print_benchmark=False, benchmark_name='') -> Tuple[np.ndarray, n self.alfa) i = 0 while self.t[i] <= self.tmax: - # 2do: correct alfa coefficients in case of changing dt according to p.38 of [Everhart1] u_next, du_dt_next, self.alfa, real_tolerance = self._do_step(self.u, self.du_dt, self.t, self.f2, @@ -919,8 +919,9 @@ def solve(self, print_benchmark=False, benchmark_name='') -> Tuple[np.ndarray, n if self.is_adaptive_step: # according to chapter 3.4 from [Everhart1] - # 0.5 is my own empirical coefficient - self.dt = 0.5 * ((self.tolerance / real_tolerance) ** (1 / ((len(self.h) - 1) + 2))) + # coeff is my own empirical coefficient + coeff = 0.5 + self.dt = coeff * ((self.tolerance / real_tolerance) ** (1 / ((len(self.h) - 1) + 2))) self.u = np.vstack([self.u, u_next]) self.du_dt = np.vstack([self.du_dt, du_dt_next]) @@ -947,9 +948,9 @@ def _do_step(self, u, du_dt, t, f, dt, h, alfa) -> Tuple[np.ndarray, np.ndarray, a = np.zeros([a_size, u_size], dtype='longdouble') # calculate c coefficients according to (9) from [Everhart1] - c = np.zeros([tau_size, tau_size], dtype='longdouble') - for i in range(tau_size): - for j in range(tau_size): + c = np.zeros([a_size, a_size], dtype='longdouble') + for i in range(a_size): + for j in range(a_size): if i == j: c[i, j] = 1 elif (j == 0) and (i > 0): @@ -970,13 +971,13 @@ def _do_step(self, u, du_dt, t, f, dt, h, alfa) -> Tuple[np.ndarray, np.ndarray, f_tau[i] = f(u_tau[i], du_dt_tau[i], t[-1] + tau[i]) # correct alfa coefficients according to (7) from [Everhart1] - for j in range(i): + for j in range(a_size): alfa[j] = self._divided_difference(j + 1, f_tau, tau) # correct a coefficients according to (8) from [Everhart1] - for j in range(i): + for j in range(a_size): a[j] = alfa[j] - for k in range(j + 1, tau_size): + for k in range(j + 1, a_size): a[j] += c[k, j] * alfa[k] # correct values of the function and its derivative at the end of dt interval diff --git a/setup.py b/setup.py index 09863a8..cdbea5d 100644 --- a/setup.py +++ b/setup.py @@ -7,7 +7,7 @@ setup( name="cleanode", - version="0.1.5", + version="0.1.6", author="Vladimir Kirievskiy", author_email="vlakir1234@gmail.com", description="Example using an embedded solver", diff --git a/system_ode_example.py b/system_ode_example.py index abc3a91..f1df1bb 100644 --- a/system_ode_example.py +++ b/system_ode_example.py @@ -92,6 +92,7 @@ def exact_f(t): dt0 = np.longdouble(0.01) tolerance = 1e-4 + is_adaptive_step = True # initial conditions: x0 = np.longdouble(0) @@ -102,7 +103,7 @@ def exact_f(t): vz0 = np.longdouble(0) u0 = np.array([x0, vx0, y0, vy0, z0, vz0], dtype='longdouble') - solver = Fehlberg45Solver(f, u0, t0, tmax, dt0, is_adaptive_step=True, tolerance=tolerance) + solver = Fehlberg45Solver(f, u0, t0, tmax, dt0, is_adaptive_step=is_adaptive_step, tolerance=tolerance) solution, time_points = solver.solve(print_benchmark=True, benchmark_name=solver.name) x_solution = solution[:, 0] y_solution = solution[:, 2] @@ -116,7 +117,7 @@ def exact_f(t): u0 = np.array([x0, y0, z0], dtype='longdouble') du_dt0 = np.array([vx0, vy0, vz0], dtype='longdouble') - solver1 = EverhartIIRadau7ODESolver(f2, u0, du_dt0, t0, tmax, dt0, is_adaptive_step=True, + solver1 = EverhartIIRadau7ODESolver(f2, u0, du_dt0, t0, tmax, dt0, is_adaptive_step=is_adaptive_step, tolerance=tolerance) solution1, d_solution1, time_points1 = solver1.solve(print_benchmark=True, benchmark_name=solver1.name) x_solution1 = solution1[:, 0] @@ -126,7 +127,12 @@ def exact_f(t): ax.plot(time_points1, x_solution1, label=solver1.name) points_number = int((tmax - t0) / dt0) - time_exact = np.linspace(t0, t0 + dt0 * points_number, (points_number + 1)) + + if is_adaptive_step: + time_exact = np.linspace(t0, t0 + dt0 * points_number, (points_number + 1)) + else: + time_exact = np.linspace(t0, t0 + dt0 * points_number, points_number + 2) + x_exact, y_exact = exact_f(time_exact) plt.plot(time_exact, x_exact, label='Exact analytical solution') @@ -145,6 +151,3 @@ def exact_f(t): figManager.window.showMaximized() plt.show() - - -