|
1 |
| -from typing import Callable |
| 1 | +"""This module defines the interface for solving initial-value problems |
| 2 | +for ordinary differential equations: |
| 3 | +
|
| 4 | + .. math:: |
| 5 | + \\frac{dy}{dt} = f(t, y), \\quad y(t_0) = y_0. |
| 6 | +
|
| 7 | +""" |
| 8 | + |
| 9 | +from collections.abc import Callable |
| 10 | +from typing import TypeAlias |
2 | 11 |
|
3 | 12 | import numpy as np
|
4 | 13 | from oif.core import (
|
|
13 | 22 | unload_impl,
|
14 | 23 | )
|
15 | 24 |
|
16 |
| -rhs_fn_t = Callable[[float, np.ndarray, np.ndarray, object], int] |
| 25 | +RhsFn: TypeAlias = Callable[[float, np.ndarray, np.ndarray, object], int] |
| 26 | +"""Signature of the right-hand side (RHS) function :math:`f(t, y)`. |
| 27 | +
|
| 28 | +The function accepts four arguments: |
| 29 | + - `t`: current time, |
| 30 | + - `y`: state vector at time :math:`t`, |
| 31 | + - `ydot`: output array to which the result of function evalutation is stored, |
| 32 | + - `user_data`: additional context (user-defined data) that |
| 33 | + must be passed to the function (e.g., parameters of the system). |
| 34 | +
|
| 35 | +""" |
17 | 36 |
|
18 | 37 |
|
19 | 38 | class IVP:
|
| 39 | + """Interface for solving initial value problems. |
| 40 | +
|
| 41 | + This class serves as a gateway to the implementations of the |
| 42 | + solvers for initial-value problems for ordinary differential equations. |
| 43 | +
|
| 44 | + Parameters |
| 45 | + ---------- |
| 46 | + impl : str |
| 47 | + Name of the desired implementation. |
| 48 | +
|
| 49 | + Examples |
| 50 | + -------- |
| 51 | +
|
| 52 | + Let's solve the following initial value problem: |
| 53 | +
|
| 54 | + .. math:: |
| 55 | + y'(t) = -y(t), \\quad y(0) = 1. |
| 56 | +
|
| 57 | + First, import the necessary modules: |
| 58 | + >>> import numpy as np |
| 59 | + >>> from oif.interfaces.ivp import IVP |
| 60 | +
|
| 61 | + Define the right-hand side function: |
| 62 | +
|
| 63 | + >>> def rhs(t, y, ydot, user_data): |
| 64 | + ... ydot[0] = -y[0] |
| 65 | + ... return 0 # No errors, optional |
| 66 | +
|
| 67 | + Now define the initial condition: |
| 68 | +
|
| 69 | + >>> y0, t0 = np.array([1.0]), 0.0 |
| 70 | +
|
| 71 | + Create an instance of the IVP solver using the implementation "jl_diffeq", |
| 72 | + which is an adapter to the `OrdinaryDiffeq.jl` Julia package: |
| 73 | +
|
| 74 | + >>> s = IVP("jl_diffeq") |
| 75 | +
|
| 76 | + We set the initial value, the right-hand side function, and the tolerance: |
| 77 | +
|
| 78 | + >>> s.set_initial_value(y0, t0) |
| 79 | + >>> s.set_rhs_fn(rhs) |
| 80 | + >>> s.set_tolerances(1e-6, 1e-12) |
| 81 | +
|
| 82 | + Now we integrate to time `t = 1.0` in a loop, outputting the current value |
| 83 | + of `y` with time step `0.1`: |
| 84 | +
|
| 85 | + >>> t = t0 |
| 86 | + >>> times = np.linspace(t0, t0 + 1.0, num=11) |
| 87 | + >>> for t in times[1:]: |
| 88 | + ... s.integrate(t) |
| 89 | + ... print(f"{t:.1f} {s.y[0]:.6f}") |
| 90 | + 0.1 0.904837 |
| 91 | + 0.2 0.818731 |
| 92 | + 0.3 0.740818 |
| 93 | + 0.4 0.670320 |
| 94 | + 0.5 0.606531 |
| 95 | + 0.6 0.548812 |
| 96 | + 0.7 0.496585 |
| 97 | + 0.8 0.449329 |
| 98 | + 0.9 0.406570 |
| 99 | + 1.0 0.367879 |
| 100 | +
|
| 101 | + """ |
| 102 | + |
20 | 103 | def __init__(self, impl: str):
|
21 | 104 | self._binding: OIFPyBinding = load_impl("ivp", impl, 1, 0)
|
22 |
| - self.s = None |
23 |
| - self.N: int = 0 |
24 |
| - self.y0: np.ndarray |
25 | 105 | self.y: np.ndarray
|
| 106 | + """Current value of the state vector.""" |
| 107 | + self._N: int = 0 |
26 | 108 |
|
27 | 109 | def set_initial_value(self, y0: np.ndarray, t0: float):
|
| 110 | + """Set initial value y(t0) = y0.""" |
28 | 111 | y0 = np.asarray(y0, dtype=np.float64)
|
29 | 112 | self.y0 = y0
|
30 | 113 | self.y = np.empty_like(y0)
|
31 |
| - self.N = len(self.y0) |
| 114 | + self._N = len(self.y0) |
32 | 115 | t0 = float(t0)
|
33 | 116 | self._binding.call("set_initial_value", (y0, t0), ())
|
34 | 117 |
|
35 |
| - def set_rhs_fn(self, rhs_fn: rhs_fn_t): |
36 |
| - if self.N <= 0: |
| 118 | + def set_rhs_fn(self, rhs_fn: RhsFn): |
| 119 | + """Specify right-hand side function f.""" |
| 120 | + if self._N <= 0: |
37 | 121 | raise RuntimeError("'set_initial_value' must be called before 'set_rhs_fn'")
|
38 | 122 |
|
39 | 123 | self.wrapper = make_oif_callback(
|
40 | 124 | rhs_fn, (OIF_FLOAT64, OIF_ARRAY_F64, OIF_ARRAY_F64, OIF_USER_DATA), OIF_INT
|
41 | 125 | )
|
42 | 126 | self._binding.call("set_rhs_fn", (self.wrapper,), ())
|
43 | 127 |
|
| 128 | + def set_tolerances(self, rtol: float, atol: float): |
| 129 | + """Specify relative and absolute tolerances, respectively.""" |
| 130 | + self._binding.call("set_tolerances", (rtol, atol), ()) |
| 131 | + |
44 | 132 | def set_user_data(self, user_data: object):
|
| 133 | + """Specify additional data that will be used for right-hand side function.""" |
45 | 134 | self.user_data = make_oif_user_data(user_data)
|
46 | 135 | self._binding.call("set_user_data", (self.user_data,), ())
|
47 | 136 |
|
48 |
| - def set_tolerances(self, rtol: float, atol: float): |
49 |
| - self._binding.call("set_tolerances", (rtol, atol), ()) |
| 137 | + def set_integrator(self, integrator_name: str, integrator_params: dict = {}): |
| 138 | + """Set integrator, if the name is recognizable.""" |
| 139 | + self._binding.call("set_integrator", (integrator_name, integrator_params), ()) |
50 | 140 |
|
51 | 141 | def integrate(self, t: float):
|
| 142 | + """Integrate to time `t` and write solution to `y`.""" |
52 | 143 | self._binding.call("integrate", (t,), (self.y,))
|
53 | 144 |
|
54 |
| - def set_integrator(self, integrator_name: str, integrator_params: dict = {}): |
55 |
| - self._binding.call("set_integrator", (integrator_name, integrator_params), ()) |
56 |
| - |
57 | 145 | def print_stats(self):
|
| 146 | + """Print integration statistics.""" |
58 | 147 | self._binding.call("print_stats", (), ())
|
59 | 148 |
|
60 | 149 | def __del__(self):
|
|
0 commit comments