v0.0.8 (#14)

Author: ngam

.devcontainer/devcontainer.json

@@ -40,7 +40,6 @@
 		"ms-vscode.vscode-github-issue-notebooks",
 		"ms-vscode.vscode-markdown-notebook",
 		"ms-python.gather",
-		"github.copilot",
 		"github.vscode-pull-request-github",
 		"visualstudioexptteam.vscodeintellicode",
 		"ms-vscode.github-issues-prs",

.github/workflows/test.yaml

@@ -14,7 +14,7 @@ jobs:
       fail-fast: false
       matrix:
          python-version: ["3.7", "3.8", "3.9", "3.10"]
-         backend: ["jax", "numpy", "numba"]
+         backend: ["numpy", "numba"]
          os: ["ubuntu-latest"]
 
     steps:

hypersolver/__init__.py

@@ -37,21 +37,27 @@
 
 import os
 
-from hypersolver.util import jxt as jit
+from hypersolver.util import xnp
+from hypersolver.util import jxt
+
 from hypersolver.lax_friedrichs import lx_loop
+
 from hypersolver.lax_wendroff import lw_loop
+
 from hypersolver.runge_kutta import rk_loop
 
 
-__version__ = "0.0.7"
+__version__ = "0.0.8"
 
 __hyper_methods__ = [
     "lax_friedrichs",
     "lax_wendroff",
-    "method_of_characteristics",
     "runge_kutta_2",
 ]
 
+np = xnp
+jit = jxt
+
 
 def set_solver(
     method=os.environ.get("HS_METHOD", "lax_friedrichs"),

hypersolver/derivative.py

@@ -2,7 +2,9 @@
 
     available functions:
         - ord1_acc2: order=1, accuracy=2
+        - ord1_acc4: order=1, accuracy=4
         - ord2_acc2: order=2, accuracy=2
+        - ord2_acc4: order=2, accuracy=4
 """
 
 from hypersolver.util import xnp as np
@@ -79,3 +81,39 @@ def ord2_acc2(_func, _xvar):
     )
 
     return np.concatenate((_derivative0, _derivative, _derivative1))
+
+
+@jit(nopython=True)
+def ord2_acc4(_func, _xvar):
+    """ central differencing: order=2, accuracy=4 """
+
+    _func, _xvar = np.asarray(_func), np.asarray(_xvar)
+
+    _result10 = (
+        _func[:-2] - 2.0*_func[1:-1] + _func[2:]
+    )/((_xvar[2:] - _xvar[:-2])/2.0)**2
+
+    _result20 = (
+        (4/3)*_func[:-2] - (5/2)*_func[1:-1] + (4/3)*_func[2:]
+    )/((_xvar[2:] - _xvar[:-2])/2.0)**2
+
+    _result2 = np.concatenate((
+        np.array([_result10[0]]),
+        (
+            -(1/12)*_func[:-4] - (1/12)*_func[4:]
+        )/((_xvar[4:] - _xvar[:-4])/4.0)**2 + _result20[1:-1],
+        np.array([_result10[-1]]),
+    ))
+
+    _derivative = _result2
+
+    _derivative0 = np.asarray(
+        [(_func[2] - 2.0*_func[1] + _func[0])/(_xvar[1] - _xvar[0])**2]
+    )
+    _derivative1 = np.asarray(
+        [(_func[-1] - 2.0*_func[-2] + _func[-3])/(_xvar[-1] - _xvar[-2])**2]
+    )
+
+    return np.concatenate(
+        (_derivative0, _derivative, _derivative1)
+    )

hypersolver/lax_wendroff.py

@@ -35,6 +35,7 @@
 from hypersolver.util import jxt as jit
 from hypersolver.util import xnp as np
 from hypersolver.util import term_util, time_step_util
+
 from hypersolver.derivative import ord1_acc2, ord2_acc2
 
 

hypersolver/runge_kutta.py

@@ -7,6 +7,7 @@
     init_: n
     vars_: x
     func_: f
+    sink_: 0
     time_: Δt
 
     outputs:
@@ -15,7 +16,7 @@
 
     numerics:
     ---------
-    - next_vals: n + Δt*f(n + Δt/2*f(n, x))
+    - rk2 next_: n + Δt*f(n + Δt/2*f(n, x))
 
 """
 

hypersolver/tests/derivative_test.py

@@ -5,7 +5,7 @@
 
 from hypersolver.util import xnp as np
 from hypersolver.derivative import ord1_acc2, ord2_acc2
-from hypersolver.derivative import ord1_acc4
+from hypersolver.derivative import ord1_acc4, ord2_acc4
 
 xvar = np.linspace(0, 10, 1000)
 yvar = np.sin(xvar)
@@ -75,3 +75,31 @@ def test_ord2_acc2():
         assert (ord2_acc2(
             yvar, xvar) - d2ydx2).sum() == pytest.approx(
                 (d2ydx2-d2ydx2).sum(), abs=1e-1)
+
+
+def test_ord2_acc4():
+    """ test: central differncing: order=2, accuracy=2 """
+
+    assert ord2_acc4(
+        yvar, xvar).shape == d2ydx2.shape
+
+    assert ord2_acc4(
+        yvar, xvar) == pytest.approx(d2ydx2, abs=1e-1)
+
+    if os.environ.get("HS_BACKEND", "numpy") == "numpy":
+        assert (ord2_acc4(
+            yvar, xvar) - d2ydx2).sum() == pytest.approx(0.0, abs=1e-1)
+    else:
+        assert (ord2_acc4(
+            yvar, xvar) - d2ydx2).sum() == pytest.approx(
+                (d2ydx2-d2ydx2).sum(), abs=1e-1)
+
+
+def test_comp_ord2():
+    """ test: compare ord1 acc(2,4) """
+
+    assert (
+        (ord2_acc4(yvar, xvar)-d2ydx2)**2
+    ).sum() <= (
+        (ord2_acc2(yvar, xvar)-d2ydx2)**2
+    ).sum()

hypersolver/util.py

@@ -14,9 +14,8 @@ def set_xnp(backend=os.environ.get("HS_BACKEND", "numpy")):
 
     if backend == "jax":
         warnings.warn(
-            "experimental jax support is suboptimal with no performance gain")
-        import jax.numpy as jnp  # pylint: disable=import-outside-toplevel
-        return jnp
+            f"no more {backend} support, reverting to numpy")
+        return np
 
     return np
 
@@ -66,7 +65,7 @@ def term_util(term, orig):
         xnp.asarray(term, dtype=orig.dtype), orig)[0]
 
 
-# @jxt(nopython=True, parallel=True)
+# @jxt(nopython=True)
 def func_util(func, _vals, _vars, **kwargs):
     """ evaluate function if one """
     return func(_vals, _vars, **kwargs) if callable(func) else func

pyproject.toml

@@ -31,7 +31,7 @@ repository = "https://github.com/ngam/hypersolver"
 dev = [
     "pylint", "pytest", "autopep8", "jupyterlab", "pytype",
     "typing", "build", "flake8", "jupyter-book", "ghp-import",
-    "matplotlib", "jax", "jaxlib", "numba", "pandas", "sympy",
+    "matplotlib", "numba",
 ]
 
 [build-system]