bugfix

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')
 

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

setup.py

@@ -7,7 +7,7 @@
 
 setup(
     name="cleanode",
-    version="0.1.5",
+    version="0.1.6",
     author="Vladimir Kirievskiy",
     author_email="[email protected]",
     description="Example using an embedded solver",

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()
-
-
-