68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324 | def solve_VIDE(*, kernel_values, a_values=None, g_values=None, soln_init_value, time_step=1.0,
coll_divs=2, coll_choices=[0,1,2], return_polys=False, show_warnings=True):
r'''
Solve a Volterra integro-differential equation.
Finds $y(t)$ satisfying
$$y'(t) = a(t)\,y(t) + g(t) + \int_0^t K(t-s)\,y(s)\,ds, \quad y(0) = y_0$$
Parameters
----------
kernel_values : array_like of shape (N,) or (N, d, d)
Values of $K(s)$ at times $s = 0, h, 2h, \ldots, (N-1)h$, where $h$
is ``time_step``. Pass a 1-D array for scalar equations or a 3-D array
of shape ``(N, d, d)`` for $d$-dimensional vector equations.
a_values : array_like of shape (N,) or (N, d, d), optional
Values of the coefficient $a(t)$ at the same times as
``kernel_values``. For vector equations $a(t)$ is a $d \times d$
matrix. Defaults to zero.
g_values : array_like of shape (N,) or (N, d) or (N, d, m), optional
Forcing term $g(t)$ sampled at the same times as ``kernel_values``.
For matrix-valued equations pass shape ``(N, d, m)`` to solve $m$
right-hand sides simultaneously. Defaults to zero.
soln_init_value : float or array_like of shape (d,) or (d, m)
Initial value $y(0) = y_0$. Required.
time_step : float, optional
Spacing $h$ between consecutive sample times. Must be positive.
Default is 1.0.
coll_divs : int, optional
Number of collocation sub-intervals per mesh interval. Must be a
positive integer. Default is 2.
coll_choices : list of int, optional
Indices selecting the collocation nodes within each sub-interval.
Each entry $k$ corresponds to the node $k / c$ where $c$ =
``coll_divs``, placed in $[0, 1]$. Entries must be distinct integers
in $\{0, 1, \ldots, \text{coll\_divs}\}$. Default is ``[0, 1, 2]``.
return_polys : bool, optional
If ``True``, also return the piecewise polynomial solution.
Default is ``False``.
show_warnings : bool, optional
If ``True`` (default), print a warning when ``kernel_values`` is
truncated or when the Numba fallback is used.
Returns
-------
soln_values : ndarray of shape (N,) or (N, d) or (N, d, m)
Solution values $y(t)$ at the same times as the input arrays.
Returned when ``return_polys=False`` (default).
(soln_values, polys) : tuple
Returned when ``return_polys=True``. ``soln_values`` is as above.
For scalar equations, ``polys`` is a list of
`numpy.polynomial.Polynomial` objects, one per mesh interval, each
mapping $t$ to the polynomial approximation of $y(t)$ on that
interval. For vector equations, each element of ``polys`` is an object
array of shape ``(d,)`` (or ``(d, m)`` for matrix equations) containing
one polynomial per component.
Notes
-----
The length $N$ of the input arrays must satisfy
$N \equiv 1 \pmod{\text{coll\_divs}^2}$. If a longer array is supplied it
is truncated to the largest conforming length and a warning is printed
(unless ``show_warnings=False``).
The solver dispatches at runtime to a D-extension routine specialised for
the given collocation setting. For scalar equations, settings not compiled
into the extension fall back to a Numba-JIT implementation (requires the
``numba`` optional dependency); a warning is printed when the fallback is
used. For vector equations only the compiled settings are supported.
References
----------
.. [1] Brunner, H. *Collocation Methods for Volterra Integral and Related
Functional Differential Equations.* Cambridge University Press, 2004.
Chapter 3, pp. 160–167.
'''
# ------------------------------------------------------------------ complex dispatch
if _cplx.is_complex(kernel_values, a_values, g_values, soln_init_value):
K_arr = np.asarray(kernel_values)
is_scalar = (K_arr.ndim == 1)
d_orig = 0 if is_scalar else K_arr.shape[1]
K_real = _cplx._block_kernel(K_arr)
a_real = _cplx._block_a(np.asarray(a_values)) if a_values is not None else None
g_real = _cplx._expand_g(np.asarray(g_values)) if g_values is not None else None
init_real = _cplx._expand_init(soln_init_value)
result = solve_VIDE(
kernel_values=K_real, a_values=a_real, g_values=g_real,
soln_init_value=init_real, time_step=time_step, coll_divs=coll_divs,
coll_choices=coll_choices, return_polys=return_polys,
show_warnings=show_warnings)
if return_polys:
soln_real, polys_real = result
return (_cplx._recombine(soln_real, d_orig),
_cplx._recombine_polys(polys_real, d_orig))
return _cplx._recombine(result, d_orig)
kernel_values_ = np.asarray(kernel_values, dtype=float)
ndim = kernel_values_.ndim
if ndim not in (1, 3):
raise ValueError(
f"kernel_values must be 1-D (scalar) or 3-D (N, d, d), got shape {kernel_values_.shape}")
N_orig = len(kernel_values_)
N, kernel_values_ = _truncate_N(kernel_values_, coll_divs, show_warnings)
# ------------------------------------------------------------------ vector path
if ndim == 3:
_, d1, d2 = kernel_values_.shape
if d1 != d2:
raise ValueError(f"kernel_values must have shape (N, d, d), got {kernel_values_.shape}")
d = d1
# ---- matrix case: detect via soln_init_value shape (d, m_cols) ----
soln_init_values_ = np.asarray(soln_init_value, dtype=float)
if soln_init_values_.ndim == 2:
d_init, m_cols = soln_init_values_.shape
if d_init != d:
raise ValueError(
f"soln_init_value shape {soln_init_values_.shape} incompatible with d={d}")
if g_values is not None:
g_mat = np.asarray(g_values, dtype=float)
if g_mat.ndim == 3:
if g_mat.shape[1:] != (d, m_cols):
raise ValueError(
f"g_values shape {g_mat.shape} incompatible with kernel/soln_init shapes")
g_cols = [g_mat[:N, :, j] for j in range(m_cols)]
else:
g_cols = [g_values] * m_cols
else:
g_cols = [None] * m_cols
a_trunc = np.asarray(a_values, dtype=float)[:N] if a_values is not None else None
def _col_vide(j):
return solve_VIDE(kernel_values=kernel_values_,
a_values=a_trunc,
g_values=g_cols[j],
soln_init_value=soln_init_values_[:, j],
time_step=time_step, coll_divs=coll_divs,
coll_choices=coll_choices,
return_polys=return_polys,
show_warnings=show_warnings)
with ThreadPoolExecutor(max_workers=m_cols) as ex:
results = list(ex.map(_col_vide, range(m_cols)))
if return_polys:
col_solns = [r[0] for r in results]
col_polys = [r[1] for r in results]
soln = np.stack(col_solns, axis=2)
mesh_divs = len(col_polys[0])
mat_polys = []
for n in range(mesh_divs):
arr = np.empty((d, m_cols), dtype=object)
for j in range(m_cols):
arr[:, j] = col_polys[j][n]
mat_polys.append(arr)
return (soln, mat_polys)
return np.stack(results, axis=2)
if g_values is not None:
g_values_ = np.asarray(g_values, dtype=float)
if g_values_.shape != (N_orig, d):
raise ValueError(
f"g_values shape {g_values_.shape} incompatible with kernel_values shape {kernel_values_.shape}")
g_values_ = g_values_[:N]
else:
g_values_ = np.zeros((N, d), dtype=float)
if a_values is not None:
a_values_ = np.asarray(a_values, dtype=float)
if a_values_.shape != (N_orig, d, d):
raise ValueError(
f"a_values shape {a_values_.shape} incompatible with kernel_values shape {kernel_values_.shape}")
a_values_ = a_values_[:N]
else:
a_values_ = np.zeros((N, d, d), dtype=float)
soln_init_values_ = soln_init_values_.ravel()
if soln_init_values_.shape != (d,):
raise ValueError(
f"soln_init_value must be a scalar or length-{d} array for d={d}")
assert coll_divs > 0, "coll_divs must be a positive integer"
assert all(isinstance(c, int) for c in coll_choices), "coll_choices must be a list of integers"
assert all(coll_choices.count(c) <= 1 for c in coll_choices), \
"all integers in coll_choices must be distinct"
for choice in coll_choices:
assert 0 <= choice <= coll_divs, "coll_choices must contain only integers from 0 to coll_divs"
coll_choices = sorted(coll_choices)
if (coll_divs, coll_choices) not in _fast_settings_VIDE:
raise RuntimeError(
f"Collocation setting (coll_divs={coll_divs}, coll_choices={coll_choices}) "
f"not supported by D extension.")
k_c = np.ascontiguousarray(kernel_values_, dtype=np.float64)
g_c = np.ascontiguousarray(g_values_, dtype=np.float64)
a_c = np.ascontiguousarray(a_values_, dtype=np.float64)
N_used = len(k_c)
mesh_divs = (N_used - 1) // coll_divs**2
soln_vals, poly_coefs = _dlang_module.solve_vide_vec_d(
g_c, k_c, a_c, soln_init_values_, time_step, coll_divs, coll_choices, return_polys)
if return_polys:
return (soln_vals, _build_vec_polys(poly_coefs, mesh_divs, coll_divs, time_step))
return soln_vals
# ------------------------------------------------------------------ scalar path
assert len(kernel_values_.shape) == 1, "kernel_values must be a 1-dim array"
if g_values is not None:
g_values_ = np.asarray(g_values, dtype=float)
assert len(g_values_.shape) == 1, "g_values must be a 1-dim array"
assert len(g_values_) == N_orig, "kernel_values and g_values must have the same length"
g_values_ = g_values_[:N]
else:
g_values_ = np.zeros(N)
if a_values is not None:
a_values_ = np.asarray(a_values, dtype=float)
assert len(a_values_.shape) == 1, "a_values must be a 1-dim array"
assert len(a_values_) == N_orig, "kernel_values and a_values must have the same length"
a_values_ = a_values_[:N]
else:
a_values_ = np.zeros(N)
assert coll_divs > 0, "coll_divs must be a positive integer"
assert all([isinstance(choice, int) for choice in coll_choices]), \
"coll_choices must be a list of integers"
assert all([coll_choices.count(c) <= 1 for c in coll_choices]), \
"all integers in coll_choices must be distinct"
for choice in coll_choices:
assert 0 <= choice <= coll_divs, "coll_choices must contain only integers from 0 to coll_divs"
coll_choices = sorted(coll_choices)
if (coll_divs, coll_choices) in _fast_settings_VIDE:
soln_vals, poly_coefs = _dlang_module.solve_vide_d(
g_values_, kernel_values_, a_values_, soln_init_value,
time_step, coll_divs, coll_choices, return_polys)
elif _numba_available:
if show_warnings:
print("warning: falling back to slower python/numba code")
soln_vals, poly_coefs = _numba_solvers.solve_VIDE_jit(
g_values_, kernel_values_, a_values_, soln_init_value,
time_step, coll_divs, coll_choices, return_polys)
else:
raise NotImplementedError(
f"Collocation setting (coll_divs={coll_divs}, coll_choices={coll_choices}) is not "
f"supported by the D extension. Install numba to enable the fallback solver, or "
f"use a supported setting (see fast_coll_settings_VIDE)."
)
if return_polys:
polys = []
for i, coefs in enumerate(poly_coefs):
domain = (i * coll_divs**2 * time_step, (i+1) * coll_divs**2 * time_step)
poly = np.polynomial.Polynomial(coefs, domain=domain, window=(0.0, 1.0), symbol='t')
poly = poly.convert(domain=domain, window=domain)
polys.append(poly)
return (soln_vals, polys)
else:
return soln_vals
|