mirror of
https://github.com/meteoinfo/MeteoInfo.git
synced 2025-12-08 20:36:05 +00:00
add solve_ivp ODE function
This commit is contained in:
wyq
parent
a3aefa95d3
commit
7594c6c1f6
1
.idea/compiler.xml
generated
1
.idea/compiler.xml
generated
@ -24,6 +24,7 @@
|
||||
<module name="meteoinfo-ui" />
|
||||
<module name="meteoinfo-math" />
|
||||
<module name="meteoinfo-common" />
|
||||
<module name="meteoinfo-jython" />
|
||||
</profile>
|
||||
</annotationProcessing>
|
||||
<bytecodeTargetLevel target="11" />
|
||||
|
||||
2
.idea/encodings.xml
generated
2
.idea/encodings.xml
generated
@ -17,6 +17,8 @@
|
||||
<file url="file://$PROJECT_DIR$/meteoinfo-geometry/src/main/resources" charset="UTF-8" />
|
||||
<file url="file://$PROJECT_DIR$/meteoinfo-image/src/main/java" charset="UTF-8" />
|
||||
<file url="file://$PROJECT_DIR$/meteoinfo-image/src/main/resources" charset="UTF-8" />
|
||||
<file url="file://$PROJECT_DIR$/meteoinfo-jython/src/main/java" charset="UTF-8" />
|
||||
<file url="file://$PROJECT_DIR$/meteoinfo-jython/src/main/resources" charset="UTF-8" />
|
||||
<file url="file://$PROJECT_DIR$/meteoinfo-lab/src/main/java" charset="UTF-8" />
|
||||
<file url="file://$PROJECT_DIR$/meteoinfo-lab/src/main/resources" charset="UTF-8" />
|
||||
<file url="file://$PROJECT_DIR$/meteoinfo-map/src/main/java" charset="UTF-8" />
|
||||
|
||||
@ -1544,7 +1544,11 @@ public class Plot2D extends AbstractPlot2D {
|
||||
ShapeTypes stype = ShapeTypes.POLYLINE;
|
||||
ls = new LegendScheme(stype);
|
||||
for (Graphic g : this.graphics.getGraphics()) {
|
||||
ls.getLegendBreaks().add(g.getLegend());
|
||||
if (g instanceof GraphicCollection) {
|
||||
ls.addLegendBreak(((GraphicCollection) g).getLegends());
|
||||
} else {
|
||||
ls.addLegendBreak(g.getLegend());
|
||||
}
|
||||
/*if (g.getShapeType() == ShapeTypes.POLYLINE) {
|
||||
ls.getLegendBreaks().add(g.getLegend());
|
||||
}*/
|
||||
|
||||
@ -69,7 +69,7 @@ import java.util.zip.ZipInputStream;
|
||||
public static String getVersion() {
|
||||
String version = GlobalUtil.class.getPackage().getImplementationVersion();
|
||||
if (version == null || version.equals("")) {
|
||||
version = "4.0.4";
|
||||
version = "4.0.5";
|
||||
}
|
||||
return version;
|
||||
}
|
||||
|
||||
@ -32,11 +32,6 @@
|
||||
<artifactId>autocomplete</artifactId>
|
||||
<version>3.3.1</version>
|
||||
</dependency>
|
||||
<dependency>
|
||||
<groupId>org.meteothink</groupId>
|
||||
<artifactId>meteoinfo-ndarray</artifactId>
|
||||
<version>${project.version}</version>
|
||||
</dependency>
|
||||
</dependencies>
|
||||
|
||||
<build>
|
||||
|
||||
@ -1,63 +0,0 @@
|
||||
package org.meteoinfo.console.jython;
|
||||
|
||||
import org.checkerframework.checker.units.qual.C;
|
||||
import org.meteoinfo.ndarray.Array;
|
||||
import org.meteoinfo.ndarray.Complex;
|
||||
import org.meteoinfo.ndarray.DataType;
|
||||
import org.python.core.PyComplex;
|
||||
|
||||
import java.util.List;
|
||||
|
||||
public class JythonUtil {
|
||||
|
||||
/**
|
||||
* Convert PyComplex list to ArrayComplex
|
||||
* @param data PyComplex list
|
||||
* @return ArrayComplex
|
||||
*/
|
||||
public static Array toComplexArray(PyComplex data) {
|
||||
Array a = Array.factory(DataType.COMPLEX, new int[]{1});
|
||||
Complex d = new Complex(data.real, data.imag);
|
||||
a.setComplex(0, d);
|
||||
return a;
|
||||
}
|
||||
|
||||
/**
|
||||
* Convert PyComplex list to ArrayComplex
|
||||
* @param data PyComplex list
|
||||
* @return ArrayComplex
|
||||
*/
|
||||
public static Array toComplexArray(List<Object> data) {
|
||||
if (data.get(0) instanceof List) {
|
||||
int ndim = data.size();
|
||||
int len = ((List) data.get(0)).size();
|
||||
Array a = Array.factory(DataType.COMPLEX, new int[]{ndim, len});
|
||||
PyComplex pd;
|
||||
for (int i = 0; i < ndim; i++) {
|
||||
List<Object> d = (List) data.get(i);
|
||||
for (int j = 0; j < len; j++) {
|
||||
if (d.get(j) instanceof PyComplex) {
|
||||
pd = (PyComplex) d.get(j);
|
||||
a.setComplex(i * len + j, new Complex(pd.real, pd.imag));
|
||||
} else {
|
||||
a.setComplex(i * len + j, new Complex((double) d.get(j), 0));
|
||||
}
|
||||
}
|
||||
}
|
||||
return a;
|
||||
} else {
|
||||
Array a = Array.factory(DataType.COMPLEX, new int[]{data.size()});
|
||||
PyComplex pd;
|
||||
for (int i = 0; i < data.size(); i++) {
|
||||
if (data.get(i) instanceof PyComplex) {
|
||||
pd = (PyComplex) data.get(i);
|
||||
a.setComplex(i, new Complex(pd.real, pd.imag));
|
||||
} else {
|
||||
a.setComplex(i, new Complex((double) data.get(i), 0));
|
||||
}
|
||||
}
|
||||
return a;
|
||||
}
|
||||
}
|
||||
|
||||
}
|
||||
@ -483,6 +483,19 @@ public class GraphicCollection extends Graphic implements Iterator {
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Get legend list
|
||||
* @return Legend list
|
||||
*/
|
||||
public List<ColorBreak> getLegends() {
|
||||
List<ColorBreak> breaks = new ArrayList<>();
|
||||
for (Graphic graphic : this.graphics) {
|
||||
breaks.add(graphic.getLegend());
|
||||
}
|
||||
|
||||
return breaks;
|
||||
}
|
||||
|
||||
/**
|
||||
* Select graphics by an extent
|
||||
*
|
||||
|
||||
@ -1,32 +1,32 @@
|
||||
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
|
||||
<MeteoInfo File="milconfig.xml" Type="configurefile">
|
||||
<Path OpenPath="D:\Working\MIScript\Jython\mis\io\geotiff">
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\satellite\calipso"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\burf"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\satellite"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\satellite\ascat"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\image"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\ndimage"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\micaps"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\zarr"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\plot"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types"/>
|
||||
<Path OpenPath="D:\Working\MIScript\Jython\mis\common_math\integrate">
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\plot_types\funny"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\meteo\calc"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\geotiff"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\ndimage"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\array\complex"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\array"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\dataframe"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\io\micaps"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\chart"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\chart\legend"/>
|
||||
<RecentFolder Folder="D:\Working\MIScript\Jython\mis\common_math\integrate"/>
|
||||
</Path>
|
||||
<File>
|
||||
<OpenedFiles>
|
||||
<OpenedFile File="D:\Working\MIScript\Jython\mis\satellite\ascat\ascat_1.py"/>
|
||||
<OpenedFile File="D:\Working\MIScript\Jython\mis\plot_types\funny\bird_1.py"/>
|
||||
<OpenedFile File="D:\Working\MIScript\Jython\mis\io\geotiff\geotif_KX10.py"/>
|
||||
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\odeint_lorenz.py"/>
|
||||
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\solve_ivp_celestial_motion.py"/>
|
||||
<OpenedFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\solve_ivp_Lotka-Volterra.py"/>
|
||||
</OpenedFiles>
|
||||
<RecentFiles>
|
||||
<RecentFile File="D:\Working\MIScript\Jython\mis\satellite\ascat\ascat_1.py"/>
|
||||
<RecentFile File="D:\Working\MIScript\Jython\mis\plot_types\funny\bird_1.py"/>
|
||||
<RecentFile File="D:\Working\MIScript\Jython\mis\io\geotiff\geotif_KX10.py"/>
|
||||
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\odeint_lorenz.py"/>
|
||||
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\solve_ivp_celestial_motion.py"/>
|
||||
<RecentFile File="D:\Working\MIScript\Jython\mis\common_math\integrate\solve_ivp_Lotka-Volterra.py"/>
|
||||
</RecentFiles>
|
||||
</File>
|
||||
<Font>
|
||||
@ -34,5 +34,5 @@
|
||||
</Font>
|
||||
<LookFeel DockWindowDecorated="true" LafDecorated="true" Name="FlatDarkLaf"/>
|
||||
<Figure DoubleBuffering="true"/>
|
||||
<Startup MainFormLocation="-6,0" MainFormSize="1339,808"/>
|
||||
<Startup MainFormLocation="-6,-6" MainFormSize="1292,764"/>
|
||||
</MeteoInfo>
|
||||
|
||||
@ -28,6 +28,11 @@
|
||||
<artifactId>meteoinfo-console</artifactId>
|
||||
<version>${project.version}</version>
|
||||
</dependency>
|
||||
<dependency>
|
||||
<groupId>${project.groupId}</groupId>
|
||||
<artifactId>meteoinfo-jython</artifactId>
|
||||
<version>${project.version}</version>
|
||||
</dependency>
|
||||
<dependency>
|
||||
<groupId>de.sciss</groupId>
|
||||
<artifactId>docking-frames-common</artifactId>
|
||||
|
||||
Binary file not shown.
@ -72,11 +72,11 @@ def jdate(t):
|
||||
|
||||
def jdatetime(t):
|
||||
"""
|
||||
Convert python date to joda DateTime.
|
||||
Convert python date to java DateTime.
|
||||
|
||||
:param t: Python date
|
||||
|
||||
:returns: Joda DateTime
|
||||
:returns: Java DateTime
|
||||
"""
|
||||
if isinstance(t, (list, tuple)):
|
||||
r = []
|
||||
|
||||
@ -1,3 +1,5 @@
|
||||
from ._ode import *
|
||||
from ._ivp import *
|
||||
|
||||
__all__ = _ode.__all__
|
||||
__all__ += _ivp.__all__
|
||||
|
||||
@ -0,0 +1,4 @@
|
||||
from .ivp import *
|
||||
|
||||
__all__ = []
|
||||
__all__ += ivp._all_
|
||||
166
meteoinfo-lab/pylib/mipylib/numeric/integrate/_ivp/common.py
Normal file
166
meteoinfo-lab/pylib/mipylib/numeric/integrate/_ivp/common.py
Normal file
@ -0,0 +1,166 @@
|
||||
|
||||
from java.lang import Double
|
||||
from org.meteoinfo.math.integrate import ODESolver
|
||||
|
||||
from warnings import warn
|
||||
from mipylib.numeric import core as np
|
||||
|
||||
|
||||
ESP = Double.MIN_VALUE
|
||||
|
||||
|
||||
def validate_first_step(first_step, t0, t_bound):
|
||||
"""Assert that first_step is valid and return it."""
|
||||
if first_step <= 0:
|
||||
raise ValueError("`first_step` must be positive.")
|
||||
if first_step > np.abs(t_bound - t0):
|
||||
raise ValueError("`first_step` exceeds bounds.")
|
||||
return first_step
|
||||
|
||||
|
||||
def validate_max_step(max_step):
|
||||
"""Assert that max_Step is valid and return it."""
|
||||
if max_step <= 0:
|
||||
raise ValueError("`max_step` must be positive.")
|
||||
return max_step
|
||||
|
||||
|
||||
def warn_extraneous(extraneous):
|
||||
"""Display a warning for extraneous keyword arguments.
|
||||
|
||||
The initializer of each solver class is expected to collect keyword
|
||||
arguments that it doesn't understand and warn about them. This function
|
||||
prints a warning for each key in the supplied dictionary.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
extraneous : dict
|
||||
Extraneous keyword arguments
|
||||
"""
|
||||
if extraneous:
|
||||
warn("The following arguments have no effect for a chosen solver: " \
|
||||
"{', '.join(f'`{x}`' for x in extraneous)}.",
|
||||
stacklevel=3)
|
||||
|
||||
|
||||
def validate_tol(rtol, atol, n):
|
||||
"""Validate tolerance values."""
|
||||
|
||||
if np.any(rtol < 100 * EPS):
|
||||
warn("At least one element of `rtol` is too small. " \
|
||||
"Setting `rtol = np.maximum(rtol, {})`.".format(100 * EPS),
|
||||
stacklevel=3)
|
||||
rtol = np.maximum(rtol, 100 * EPS)
|
||||
|
||||
atol = np.asarray(atol)
|
||||
if atol.ndim > 0 and atol.shape != (n,):
|
||||
raise ValueError("`atol` has wrong shape.")
|
||||
|
||||
if np.any(atol < 0):
|
||||
raise ValueError("`atol` must be positive.")
|
||||
|
||||
return rtol, atol
|
||||
|
||||
|
||||
def norm(x):
|
||||
"""Compute RMS norm."""
|
||||
return np.linalg.norm(x) / x.size ** 0.5
|
||||
|
||||
|
||||
def select_initial_step(fun, t0, y0, t_bound,
|
||||
max_step, f0, direction, order, rtol, atol):
|
||||
"""Empirically select a good initial step.
|
||||
|
||||
The algorithm is described in [1]_.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
fun : callable
|
||||
Right-hand side of the system.
|
||||
t0 : float
|
||||
Initial value of the independent variable.
|
||||
y0 : ndarray, shape (n,)
|
||||
Initial value of the dependent variable.
|
||||
t_bound : float
|
||||
End-point of integration interval; used to ensure that t0+step<=tbound
|
||||
and that fun is only evaluated in the interval [t0,tbound]
|
||||
max_step : float
|
||||
Maximum allowable step size.
|
||||
f0 : ndarray, shape (n,)
|
||||
Initial value of the derivative, i.e., ``fun(t0, y0)``.
|
||||
direction : float
|
||||
Integration direction.
|
||||
order : float
|
||||
Error estimator order. It means that the error controlled by the
|
||||
algorithm is proportional to ``step_size ** (order + 1)`.
|
||||
rtol : float
|
||||
Desired relative tolerance.
|
||||
atol : float
|
||||
Desired absolute tolerance.
|
||||
|
||||
Returns
|
||||
-------
|
||||
h_abs : float
|
||||
Absolute value of the suggested initial step.
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
|
||||
Equations I: Nonstiff Problems", Sec. II.4.
|
||||
"""
|
||||
if y0.size == 0:
|
||||
return np.inf
|
||||
|
||||
interval_length = abs(t_bound - t0)
|
||||
if interval_length == 0.0:
|
||||
return 0.0
|
||||
|
||||
scale = atol + np.abs(y0) * rtol
|
||||
d0 = norm(y0 / scale)
|
||||
d1 = norm(f0 / scale)
|
||||
if d0 < 1e-5 or d1 < 1e-5:
|
||||
h0 = 1e-6
|
||||
else:
|
||||
h0 = 0.01 * d0 / d1
|
||||
# Check t0+h0*direction doesn't take us beyond t_bound
|
||||
h0 = min(h0, interval_length)
|
||||
y1 = y0 + h0 * direction * f0
|
||||
f1 = fun(t0 + h0 * direction, y1)
|
||||
d2 = norm((f1 - f0) / scale) / h0
|
||||
|
||||
if d1 <= 1e-15 and d2 <= 1e-15:
|
||||
h1 = max(1e-6, h0 * 1e-3)
|
||||
else:
|
||||
h1 = (0.01 / max(d1, d2)) ** (1 / (order + 1))
|
||||
|
||||
return min(100 * h0, h1, interval_length, max_step)
|
||||
|
||||
|
||||
class OdeSolution():
|
||||
"""Continuous ODE solution.
|
||||
|
||||
It is organized as a collection of `DenseOutput` objects which represent
|
||||
local interpolants. It provides an algorithm to select a right interpolant
|
||||
for each given point.
|
||||
"""
|
||||
def __init__(self, solver):
|
||||
self.solver = solver
|
||||
|
||||
def __call__(self, t):
|
||||
"""Evaluate the solution.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
t : float or array_like with shape (n_points,)
|
||||
Points to evaluate at.
|
||||
|
||||
Returns
|
||||
-------
|
||||
y : ndarray, shape (n_states,) or (n_states, n_points)
|
||||
Computed values. Shape depends on whether `t` is a scalar or a
|
||||
1-D array.
|
||||
"""
|
||||
t = np.asarray(t)
|
||||
ys = self.solver.solve(t._array)
|
||||
|
||||
return np.array(ys)
|
||||
444
meteoinfo-lab/pylib/mipylib/numeric/integrate/_ivp/ivp.py
Normal file
444
meteoinfo-lab/pylib/mipylib/numeric/integrate/_ivp/ivp.py
Normal file
@ -0,0 +1,444 @@
|
||||
|
||||
from org.meteoinfo.math.integrate import ODEEquations, ODESolver
|
||||
|
||||
from mipylib.numeric import core as np
|
||||
from ...optimize import OptimizeResult
|
||||
from .common import OdeSolution
|
||||
|
||||
|
||||
_all_ = ['solve_ivp']
|
||||
|
||||
|
||||
METHODS = ['RK23', 'RK45', 'DOP853', 'Radau', 'BDF', 'LSODA']
|
||||
|
||||
|
||||
class OdeResult(OptimizeResult):
|
||||
pass
|
||||
|
||||
class ODE(ODEEquations):
|
||||
|
||||
def __init__(self, f):
|
||||
"""
|
||||
Initialize
|
||||
|
||||
:param f: Jython function
|
||||
"""
|
||||
self.f = f
|
||||
self._args = list(f.__code__.co_varnames)[2:]
|
||||
self._args = tuple(self._args)
|
||||
self.order = len(self._args)
|
||||
|
||||
def doComputeDerivatives(self, y, t):
|
||||
y = np.array(y)
|
||||
args = tuple(self.getParameters())
|
||||
r = self.f(t, y, *args)
|
||||
if isinstance(r, np.NDArray):
|
||||
r = tuple(r.aslist())
|
||||
return r
|
||||
|
||||
def solve_ivp(fun, t_span, y0, method='RK45', t_eval=None, dense_output=False, args=None, **options):
|
||||
"""Solve an initial value problem for a system of ODEs.
|
||||
|
||||
This function numerically integrates a system of ordinary differential
|
||||
equations given an initial value::
|
||||
|
||||
dy / dt = f(t, y)
|
||||
y(t0) = y0
|
||||
|
||||
Here t is a 1-D independent variable (time), y(t) is an
|
||||
N-D vector-valued function (state), and an N-D
|
||||
vector-valued function f(t, y) determines the differential equations.
|
||||
The goal is to find y(t) approximately satisfying the differential
|
||||
equations, given an initial value y(t0)=y0.
|
||||
|
||||
Some of the solvers support integration in the complex domain, but note
|
||||
that for stiff ODE solvers, the right-hand side must be
|
||||
complex-differentiable (satisfy Cauchy-Riemann equations [11]_).
|
||||
To solve a problem in the complex domain, pass y0 with a complex data type.
|
||||
Another option always available is to rewrite your problem for real and
|
||||
imaginary parts separately.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
fun : callable
|
||||
Right-hand side of the system: the time derivative of the state ``y``
|
||||
at time ``t``. The calling signature is ``fun(t, y)``, where ``t`` is a
|
||||
scalar and ``y`` is a ndarray with ``len(y) = len(y0)``. Additional
|
||||
arguments need to be passed if ``args`` is used (see documentation of
|
||||
``args`` argument). ``fun`` must return an array of the same shape as
|
||||
``y``. See `vectorized` for more information.
|
||||
t_span : 2-member sequence
|
||||
Interval of integration (t0, tf). The solver starts with t=t0 and
|
||||
integrates until it reaches t=tf. Both t0 and tf must be floats
|
||||
or values interpretable by the float conversion function.
|
||||
y0 : array_like, shape (n,)
|
||||
Initial state. For problems in the complex domain, pass `y0` with a
|
||||
complex data type (even if the initial value is purely real).
|
||||
method : string or `OdeSolver`, optional
|
||||
Integration method to use:
|
||||
|
||||
* 'RK45' (default): Explicit Runge-Kutta method of order 5(4) [1]_.
|
||||
The error is controlled assuming accuracy of the fourth-order
|
||||
method, but steps are taken using the fifth-order accurate
|
||||
formula (local extrapolation is done). A quartic interpolation
|
||||
polynomial is used for the dense output [2]_. Can be applied in
|
||||
the complex domain.
|
||||
* 'RK23': Explicit Runge-Kutta method of order 3(2) [3]_. The error
|
||||
is controlled assuming accuracy of the second-order method, but
|
||||
steps are taken using the third-order accurate formula (local
|
||||
extrapolation is done). A cubic Hermite polynomial is used for the
|
||||
dense output. Can be applied in the complex domain.
|
||||
* 'DOP853': Explicit Runge-Kutta method of order 8 [13]_.
|
||||
Python implementation of the "DOP853" algorithm originally
|
||||
written in Fortran [14]_. A 7-th order interpolation polynomial
|
||||
accurate to 7-th order is used for the dense output.
|
||||
Can be applied in the complex domain.
|
||||
* 'Radau': Implicit Runge-Kutta method of the Radau IIA family of
|
||||
order 5 [4]_. The error is controlled with a third-order accurate
|
||||
embedded formula. A cubic polynomial which satisfies the
|
||||
collocation conditions is used for the dense output.
|
||||
* 'BDF': Implicit multi-step variable-order (1 to 5) method based
|
||||
on a backward differentiation formula for the derivative
|
||||
approximation [5]_. The implementation follows the one described
|
||||
in [6]_. A quasi-constant step scheme is used and accuracy is
|
||||
enhanced using the NDF modification. Can be applied in the
|
||||
complex domain.
|
||||
* 'LSODA': Adams/BDF method with automatic stiffness detection and
|
||||
switching [7]_, [8]_. This is a wrapper of the Fortran solver
|
||||
from ODEPACK.
|
||||
|
||||
Explicit Runge-Kutta methods ('RK23', 'RK45', 'DOP853') should be used
|
||||
for non-stiff problems and implicit methods ('Radau', 'BDF') for
|
||||
stiff problems [9]_. Among Runge-Kutta methods, 'DOP853' is recommended
|
||||
for solving with high precision (low values of `rtol` and `atol`).
|
||||
|
||||
If not sure, first try to run 'RK45'. If it makes unusually many
|
||||
iterations, diverges, or fails, your problem is likely to be stiff and
|
||||
you should use 'Radau' or 'BDF'. 'LSODA' can also be a good universal
|
||||
choice, but it might be somewhat less convenient to work with as it
|
||||
wraps old Fortran code.
|
||||
|
||||
You can also pass an arbitrary class derived from `OdeSolver` which
|
||||
implements the solver.
|
||||
t_eval : array_like or None, optional
|
||||
Times at which to store the computed solution, must be sorted and lie
|
||||
within `t_span`. If None (default), use points selected by the solver.
|
||||
dense_output : bool, optional
|
||||
Whether to compute a continuous solution. Default is False.
|
||||
args : tuple, optional
|
||||
Additional arguments to pass to the user-defined functions. If given,
|
||||
the additional arguments are passed to all user-defined functions.
|
||||
So if, for example, `fun` has the signature ``fun(t, y, a, b, c)``,
|
||||
then `jac` (if given) and any event functions must have the same
|
||||
signature, and `args` must be a tuple of length 3.
|
||||
first_step : float or None, optional
|
||||
Initial step size. Default is `None` which means that the algorithm
|
||||
should choose.
|
||||
max_step : float, optional
|
||||
Maximum allowed step size. Default is np.inf, i.e., the step size is not
|
||||
bounded and determined solely by the solver.
|
||||
rtol, atol : float or array_like, optional
|
||||
Relative and absolute tolerances. The solver keeps the local error
|
||||
estimates less than ``atol + rtol * abs(y)``. Here `rtol` controls a
|
||||
relative accuracy (number of correct digits), while `atol` controls
|
||||
absolute accuracy (number of correct decimal places). To achieve the
|
||||
desired `rtol`, set `atol` to be smaller than the smallest value that
|
||||
can be expected from ``rtol * abs(y)`` so that `rtol` dominates the
|
||||
allowable error. If `atol` is larger than ``rtol * abs(y)`` the
|
||||
number of correct digits is not guaranteed. Conversely, to achieve the
|
||||
desired `atol` set `rtol` such that ``rtol * abs(y)`` is always smaller
|
||||
than `atol`. If components of y have different scales, it might be
|
||||
beneficial to set different `atol` values for different components by
|
||||
passing array_like with shape (n,) for `atol`. Default values are
|
||||
1e-3 for `rtol` and 1e-6 for `atol`.
|
||||
jac : array_like, sparse_matrix, callable or None, optional
|
||||
Jacobian matrix of the right-hand side of the system with respect
|
||||
to y, required by the 'Radau', 'BDF' and 'LSODA' method. The
|
||||
Jacobian matrix has shape (n, n) and its element (i, j) is equal to
|
||||
``d f_i / d y_j``. There are three ways to define the Jacobian:
|
||||
|
||||
* If array_like or sparse_matrix, the Jacobian is assumed to
|
||||
be constant. Not supported by 'LSODA'.
|
||||
* If callable, the Jacobian is assumed to depend on both
|
||||
t and y; it will be called as ``jac(t, y)``, as necessary.
|
||||
Additional arguments have to be passed if ``args`` is
|
||||
used (see documentation of ``args`` argument).
|
||||
For 'Radau' and 'BDF' methods, the return value might be a
|
||||
sparse matrix.
|
||||
* If None (default), the Jacobian will be approximated by
|
||||
finite differences.
|
||||
|
||||
It is generally recommended to provide the Jacobian rather than
|
||||
relying on a finite-difference approximation.
|
||||
jac_sparsity : array_like, sparse matrix or None, optional
|
||||
Defines a sparsity structure of the Jacobian matrix for a finite-
|
||||
difference approximation. Its shape must be (n, n). This argument
|
||||
is ignored if `jac` is not `None`. If the Jacobian has only few
|
||||
non-zero elements in *each* row, providing the sparsity structure
|
||||
will greatly speed up the computations [10]_. A zero entry means that
|
||||
a corresponding element in the Jacobian is always zero. If None
|
||||
(default), the Jacobian is assumed to be dense.
|
||||
Not supported by 'LSODA', see `lband` and `uband` instead.
|
||||
lband, uband : int or None, optional
|
||||
Parameters defining the bandwidth of the Jacobian for the 'LSODA'
|
||||
method, i.e., ``jac[i, j] != 0 only for i - lband <= j <= i + uband``.
|
||||
Default is None. Setting these requires your jac routine to return the
|
||||
Jacobian in the packed format: the returned array must have ``n``
|
||||
columns and ``uband + lband + 1`` rows in which Jacobian diagonals are
|
||||
written. Specifically ``jac_packed[uband + i - j , j] = jac[i, j]``.
|
||||
The same format is used in `scipy.linalg.solve_banded` (check for an
|
||||
illustration). These parameters can be also used with ``jac=None`` to
|
||||
reduce the number of Jacobian elements estimated by finite differences.
|
||||
min_step : float, optional
|
||||
The minimum allowed step size for 'LSODA' method.
|
||||
By default `min_step` is zero.
|
||||
|
||||
Returns
|
||||
-------
|
||||
Bunch object with the following fields defined:
|
||||
t : ndarray, shape (n_points,)
|
||||
Time points.
|
||||
y : ndarray, shape (n, n_points)
|
||||
Values of the solution at `t`.
|
||||
sol : `OdeSolution` or None
|
||||
Found solution as `OdeSolution` instance; None if `dense_output` was
|
||||
set to False.
|
||||
t_events : list of ndarray or None
|
||||
Contains for each event type a list of arrays at which an event of
|
||||
that type event was detected. None if `events` was None.
|
||||
y_events : list of ndarray or None
|
||||
For each value of `t_events`, the corresponding value of the solution.
|
||||
None if `events` was None.
|
||||
nfev : int
|
||||
Number of evaluations of the right-hand side.
|
||||
njev : int
|
||||
Number of evaluations of the Jacobian.
|
||||
nlu : int
|
||||
Number of LU decompositions.
|
||||
status : int
|
||||
Reason for algorithm termination:
|
||||
|
||||
* -1: Integration step failed.
|
||||
* 0: The solver successfully reached the end of `tspan`.
|
||||
* 1: A termination event occurred.
|
||||
|
||||
message : string
|
||||
Human-readable description of the termination reason.
|
||||
success : bool
|
||||
True if the solver reached the interval end or a termination event
|
||||
occurred (``status >= 0``).
|
||||
|
||||
References
|
||||
----------
|
||||
.. [1] J. R. Dormand, P. J. Prince, "A family of embedded Runge-Kutta
|
||||
formulae", Journal of Computational and Applied Mathematics, Vol. 6,
|
||||
No. 1, pp. 19-26, 1980.
|
||||
.. [2] L. W. Shampine, "Some Practical Runge-Kutta Formulas", Mathematics
|
||||
of Computation,, Vol. 46, No. 173, pp. 135-150, 1986.
|
||||
.. [3] P. Bogacki, L.F. Shampine, "A 3(2) Pair of Runge-Kutta Formulas",
|
||||
Appl. Math. Lett. Vol. 2, No. 4. pp. 321-325, 1989.
|
||||
.. [4] E. Hairer, G. Wanner, "Solving Ordinary Differential Equations II:
|
||||
Stiff and Differential-Algebraic Problems", Sec. IV.8.
|
||||
.. [5] `Backward Differentiation Formula
|
||||
<https://en.wikipedia.org/wiki/Backward_differentiation_formula>`_
|
||||
on Wikipedia.
|
||||
.. [6] L. F. Shampine, M. W. Reichelt, "THE MATLAB ODE SUITE", SIAM J. SCI.
|
||||
COMPUTE., Vol. 18, No. 1, pp. 1-22, January 1997.
|
||||
.. [7] A. C. Hindmarsh, "ODEPACK, A Systematized Collection of ODE
|
||||
Solvers," IMACS Transactions on Scientific Computation, Vol 1.,
|
||||
pp. 55-64, 1983.
|
||||
.. [8] L. Petzold, "Automatic selection of methods for solving stiff and
|
||||
nonstiff systems of ordinary differential equations", SIAM Journal
|
||||
on Scientific and Statistical Computing, Vol. 4, No. 1, pp. 136-148,
|
||||
1983.
|
||||
.. [9] `Stiff equation <https://en.wikipedia.org/wiki/Stiff_equation>`_ on
|
||||
Wikipedia.
|
||||
.. [10] A. Curtis, M. J. D. Powell, and J. Reid, "On the estimation of
|
||||
sparse Jacobian matrices", Journal of the Institute of Mathematics
|
||||
and its Applications, 13, pp. 117-120, 1974.
|
||||
.. [11] `Cauchy-Riemann equations
|
||||
<https://en.wikipedia.org/wiki/Cauchy-Riemann_equations>`_ on
|
||||
Wikipedia.
|
||||
.. [12] `Lotka-Volterra equations
|
||||
<https://en.wikipedia.org/wiki/Lotka%E2%80%93Volterra_equations>`_
|
||||
on Wikipedia.
|
||||
.. [13] E. Hairer, S. P. Norsett G. Wanner, "Solving Ordinary Differential
|
||||
Equations I: Nonstiff Problems", Sec. II.
|
||||
.. [14] `Page with original Fortran code of DOP853
|
||||
<http://www.unige.ch/~hairer/software.html>`_.
|
||||
|
||||
Examples
|
||||
--------
|
||||
Basic exponential decay showing automatically chosen time points.
|
||||
|
||||
>>> from mipylib.numeric.integrate import solve_ivp
|
||||
>>> def exponential_decay(t, y): return -0.5 * y
|
||||
>>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8])
|
||||
>>> print(sol.t)
|
||||
[ 0. 0.11487653 1.26364188 3.06061781 4.81611105 6.57445806
|
||||
8.33328988 10. ]
|
||||
>>> print(sol.y)
|
||||
[[2. 1.88836035 1.06327177 0.43319312 0.18017253 0.07483045
|
||||
0.03107158 0.01350781]
|
||||
[4. 3.7767207 2.12654355 0.86638624 0.36034507 0.14966091
|
||||
0.06214316 0.02701561]
|
||||
[8. 7.5534414 4.25308709 1.73277247 0.72069014 0.29932181
|
||||
0.12428631 0.05403123]]
|
||||
|
||||
Specifying points where the solution is desired.
|
||||
|
||||
>>> sol = solve_ivp(exponential_decay, [0, 10], [2, 4, 8],
|
||||
... t_eval=[0, 1, 2, 4, 10])
|
||||
>>> print(sol.t)
|
||||
[ 0 1 2 4 10]
|
||||
>>> print(sol.y)
|
||||
[[2. 1.21305369 0.73534021 0.27066736 0.01350938]
|
||||
[4. 2.42610739 1.47068043 0.54133472 0.02701876]
|
||||
[8. 4.85221478 2.94136085 1.08266944 0.05403753]]
|
||||
|
||||
Cannon fired upward with terminal event upon impact. The ``terminal`` and
|
||||
``direction`` fields of an event are applied by monkey patching a function.
|
||||
Here ``y[0]`` is position and ``y[1]`` is velocity. The projectile starts
|
||||
at position 0 with velocity +10. Note that the integration never reaches
|
||||
t=100 because the event is terminal.
|
||||
|
||||
>>> def upward_cannon(t, y): return [y[1], -0.5]
|
||||
>>> def hit_ground(t, y): return y[0]
|
||||
>>> hit_ground.terminal = True
|
||||
>>> hit_ground.direction = -1
|
||||
>>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10], events=hit_ground)
|
||||
>>> print(sol.t_events)
|
||||
[array([40.])]
|
||||
>>> print(sol.t)
|
||||
[0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
|
||||
1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
|
||||
|
||||
Use `dense_output` and `events` to find position, which is 100, at the apex
|
||||
of the cannonball's trajectory. Apex is not defined as terminal, so both
|
||||
apex and hit_ground are found. There is no information at t=20, so the sol
|
||||
attribute is used to evaluate the solution. The sol attribute is returned
|
||||
by setting ``dense_output=True``. Alternatively, the `y_events` attribute
|
||||
can be used to access the solution at the time of the event.
|
||||
|
||||
>>> def apex(t, y): return y[1]
|
||||
>>> sol = solve_ivp(upward_cannon, [0, 100], [0, 10],
|
||||
... events=(hit_ground, apex), dense_output=True)
|
||||
>>> print(sol.t_events)
|
||||
[array([40.]), array([20.])]
|
||||
>>> print(sol.t)
|
||||
[0.00000000e+00 9.99900010e-05 1.09989001e-03 1.10988901e-02
|
||||
1.11088891e-01 1.11098890e+00 1.11099890e+01 4.00000000e+01]
|
||||
>>> print(sol.sol(sol.t_events[1][0]))
|
||||
[100. 0.]
|
||||
>>> print(sol.y_events)
|
||||
[array([[-5.68434189e-14, -1.00000000e+01]]),
|
||||
array([[1.00000000e+02, 1.77635684e-15]])]
|
||||
|
||||
As an example of a system with additional parameters, we'll implement
|
||||
the Lotka-Volterra equations [12]_.
|
||||
|
||||
>>> def lotkavolterra(t, z, a, b, c, d):
|
||||
... x, y = z
|
||||
... return [a*x - b*x*y, -c*y + d*x*y]
|
||||
...
|
||||
|
||||
We pass in the parameter values a=1.5, b=1, c=3 and d=1 with the `args`
|
||||
argument.
|
||||
|
||||
>>> sol = solve_ivp(lotkavolterra, [0, 15], [10, 5], args=(1.5, 1, 3, 1),
|
||||
... dense_output=True)
|
||||
|
||||
Compute a dense solution and plot it.
|
||||
|
||||
>>> t = np.linspace(0, 15, 300)
|
||||
>>> z = sol.sol(t)
|
||||
>>> import matplotlib.pyplot as plt
|
||||
>>> plt.plot(t, z.T)
|
||||
>>> plt.xlabel('t')
|
||||
>>> plt.legend(['x', 'y'], shadow=True)
|
||||
>>> plt.title('Lotka-Volterra System')
|
||||
>>> plt.show()
|
||||
|
||||
A couple examples of using solve_ivp to solve the differential
|
||||
equation ``y' = Ay`` with complex matrix ``A``.
|
||||
|
||||
>>> A = np.array([[-0.25 + 0.14j, 0, 0.33 + 0.44j],
|
||||
... [0.25 + 0.58j, -0.2 + 0.14j, 0],
|
||||
... [0, 0.2 + 0.4j, -0.1 + 0.97j]])
|
||||
|
||||
Solving an IVP with ``A`` from above and ``y`` as 3x1 vector:
|
||||
|
||||
>>> def deriv_vec(t, y):
|
||||
... return A @ y
|
||||
>>> result = solve_ivp(deriv_vec, [0, 25],
|
||||
... np.array([10 + 0j, 20 + 0j, 30 + 0j]),
|
||||
... t_eval=np.linspace(0, 25, 101))
|
||||
>>> print(result.y[:, 0])
|
||||
[10.+0.j 20.+0.j 30.+0.j]
|
||||
>>> print(result.y[:, -1])
|
||||
[18.46291039+45.25653651j 10.01569306+36.23293216j
|
||||
-4.98662741+80.07360388j]
|
||||
|
||||
Solving an IVP with ``A`` from above with ``y`` as 3x3 matrix :
|
||||
|
||||
>>> def deriv_mat(t, y):
|
||||
... return (A @ y.reshape(3, 3)).flatten()
|
||||
>>> y0 = np.array([[2 + 0j, 3 + 0j, 4 + 0j],
|
||||
... [5 + 0j, 6 + 0j, 7 + 0j],
|
||||
... [9 + 0j, 34 + 0j, 78 + 0j]])
|
||||
|
||||
>>> result = solve_ivp(deriv_mat, [0, 25], y0.flatten(),
|
||||
... t_eval=np.linspace(0, 25, 101))
|
||||
>>> print(result.y[:, 0].reshape(3, 3))
|
||||
[[ 2.+0.j 3.+0.j 4.+0.j]
|
||||
[ 5.+0.j 6.+0.j 7.+0.j]
|
||||
[ 9.+0.j 34.+0.j 78.+0.j]]
|
||||
>>> print(result.y[:, -1].reshape(3, 3))
|
||||
[[ 5.67451179 +12.07938445j 17.2888073 +31.03278837j
|
||||
37.83405768 +63.25138759j]
|
||||
[ 3.39949503 +11.82123994j 21.32530996 +44.88668871j
|
||||
53.17531184+103.80400411j]
|
||||
[ -2.26105874 +22.19277664j -15.1255713 +70.19616341j
|
||||
-38.34616845+153.29039931j]]
|
||||
|
||||
|
||||
"""
|
||||
if method not in METHODS:
|
||||
raise ValueError("`method` must be one of {METHODS}.")
|
||||
|
||||
t0, tf = map(float, t_span)
|
||||
first_step = options.pop('first_step', None)
|
||||
max_step = options.pop('max_step', None)
|
||||
rtol = options.pop('rtol', None)
|
||||
atol = options.pop('atol', None)
|
||||
|
||||
fun = ODE(fun)
|
||||
if args is not None:
|
||||
fun.setParameters(args)
|
||||
|
||||
if isinstance(y0, (tuple, list)):
|
||||
y0 = np.array(y0)
|
||||
|
||||
y0 = y0.astype('double')
|
||||
ndim = len(y0)
|
||||
fun.setDimension(ndim)
|
||||
|
||||
if t_eval is not None: t_eval = np.array(t_eval)._array
|
||||
|
||||
solver = ODESolver(method, fun, t0, tf, y0._array)
|
||||
solver.setDenseOutput(dense_output)
|
||||
if t_eval is not None: solver.setTEval(t_eval)
|
||||
if first_step is not None: solver.setMinStep(first_step)
|
||||
if max_step is not None: solver.setMaxStep(max_step)
|
||||
if rtol is not None: solver.setRTol(rtol)
|
||||
if atol is not None: solver.setATol(atol)
|
||||
solver.solve()
|
||||
|
||||
ts = solver.getTResult()
|
||||
ys = solver.getYResult()
|
||||
if ts is not None: ts = np.array(ts)
|
||||
if ys is not None: ys = np.array(ys)
|
||||
|
||||
sol = OdeSolution(solver) if dense_output else None
|
||||
|
||||
return OdeResult(t=ts, y=ys, sol=sol)
|
||||
Binary file not shown.
@ -149,4 +149,48 @@ def check_random_state(seed):
|
||||
return seed
|
||||
|
||||
raise ValueError('%r cannot be used to seed a numpy.random.RandomState'
|
||||
' instance' % seed)
|
||||
' instance' % seed)
|
||||
|
||||
class _RichResult(dict):
|
||||
""" Container for multiple outputs with pretty-printing """
|
||||
def __getattr__(self, name):
|
||||
try:
|
||||
return self[name]
|
||||
except KeyError as e:
|
||||
raise AttributeError(name)
|
||||
|
||||
__setattr__ = dict.__setitem__ # type: ignore[assignment]
|
||||
__delattr__ = dict.__delitem__ # type: ignore[assignment]
|
||||
|
||||
def __repr__(self):
|
||||
order_keys = ['message', 'success', 'status', 'fun', 'funl', 'x', 'xl',
|
||||
'col_ind', 'nit', 'lower', 'upper', 'eqlin', 'ineqlin',
|
||||
'converged', 'flag', 'function_calls', 'iterations',
|
||||
'root']
|
||||
order_keys = getattr(self, '_order_keys', order_keys)
|
||||
# 'slack', 'con' are redundant with residuals
|
||||
# 'crossover_nit' is probably not interesting to most users
|
||||
omit_keys = {'slack', 'con', 'crossover_nit', '_order_keys'}
|
||||
|
||||
def key(item):
|
||||
try:
|
||||
return order_keys.index(item[0].lower())
|
||||
except ValueError: # item not in list
|
||||
return np.inf
|
||||
|
||||
def omit_redundant(items):
|
||||
for item in items:
|
||||
if item[0] in omit_keys:
|
||||
continue
|
||||
yield item
|
||||
|
||||
def item_sorter(d):
|
||||
return sorted(omit_redundant(d.items()), key=key)
|
||||
|
||||
if self.keys():
|
||||
return _dict_formatter(self, sorter=item_sorter)
|
||||
else:
|
||||
return self.__class__.__name__ + "()"
|
||||
|
||||
def __dir__(self):
|
||||
return list(self.keys())
|
||||
|
||||
Binary file not shown.
Binary file not shown.
@ -1,4 +1,5 @@
|
||||
from .minpack import *
|
||||
from ._lsq import *
|
||||
from ._optimize import *
|
||||
|
||||
__all__ = [s for s in dir() if not s.startswith('_')]
|
||||
50
meteoinfo-lab/pylib/mipylib/numeric/optimize/_optimize.py
Normal file
50
meteoinfo-lab/pylib/mipylib/numeric/optimize/_optimize.py
Normal file
@ -0,0 +1,50 @@
|
||||
|
||||
from ..lib._util import _RichResult
|
||||
|
||||
|
||||
__all__ = ['OptimizeResult']
|
||||
|
||||
|
||||
class OptimizeResult(_RichResult):
|
||||
"""
|
||||
Represents the optimization result.
|
||||
|
||||
Attributes
|
||||
----------
|
||||
x : ndarray
|
||||
The solution of the optimization.
|
||||
success : bool
|
||||
Whether or not the optimizer exited successfully.
|
||||
status : int
|
||||
Termination status of the optimizer. Its value depends on the
|
||||
underlying solver. Refer to `message` for details.
|
||||
message : str
|
||||
Description of the cause of the termination.
|
||||
fun : float
|
||||
Value of objective function at `x`.
|
||||
jac, hess : ndarray
|
||||
Values of objective function's Jacobian and its Hessian at `x` (if
|
||||
available). The Hessian may be an approximation, see the documentation
|
||||
of the function in question.
|
||||
hess_inv : object
|
||||
Inverse of the objective function's Hessian; may be an approximation.
|
||||
Not available for all solvers. The type of this attribute may be
|
||||
either np.ndarray or scipy.sparse.linalg.LinearOperator.
|
||||
nfev, njev, nhev : int
|
||||
Number of evaluations of the objective functions and of its
|
||||
Jacobian and Hessian.
|
||||
nit : int
|
||||
Number of iterations performed by the optimizer.
|
||||
maxcv : float
|
||||
The maximum constraint violation.
|
||||
|
||||
Notes
|
||||
-----
|
||||
Depending on the specific solver being used, `OptimizeResult` may
|
||||
not have all attributes listed here, and they may have additional
|
||||
attributes not listed here. Since this class is essentially a
|
||||
subclass of dict with attribute accessors, one can see which
|
||||
attributes are available using the `OptimizeResult.keys` method.
|
||||
|
||||
"""
|
||||
pass
|
||||
Binary file not shown.
@ -3836,6 +3836,7 @@ class Axes(object):
|
||||
clegend = ChartLegend(ols)
|
||||
else:
|
||||
clegend = self._axes.getLegend()
|
||||
|
||||
ls = kwargs.pop('legend', None)
|
||||
if len(args) > 0:
|
||||
if isinstance(args[0], MILayer):
|
||||
@ -3848,6 +3849,7 @@ class Axes(object):
|
||||
if not args[0].isSingleLegend():
|
||||
ls = args[0].getLegendScheme()
|
||||
args = args[1:]
|
||||
|
||||
if ls is None:
|
||||
if len(args) > 0:
|
||||
lbs = []
|
||||
@ -3858,10 +3860,12 @@ class Axes(object):
|
||||
lbs.extend(lb.legend().getLegendBreaks())
|
||||
else:
|
||||
lbs.append(lb)
|
||||
|
||||
if len(args) == 2:
|
||||
labels = args[1]
|
||||
for i in range(0, len(lbs)):
|
||||
lbs[i].setCaption(labels[i])
|
||||
|
||||
if isinstance(lbs[0], basestring):
|
||||
clegend.setTickCaptions(lbs)
|
||||
else:
|
||||
|
||||
@ -1,6 +1,6 @@
|
||||
<?xml version="1.0" encoding="UTF-8" standalone="no"?>
|
||||
<MeteoInfo File="config.xml" Type="configurefile">
|
||||
<Path OpenPath="D:\Temp\nc\2021_daily_4_VOD"/>
|
||||
<Path OpenPath="D:\Temp\image"/>
|
||||
<Font>
|
||||
<TextFont FontName="YaHei Consolas Hybrid" FontSize="14"/>
|
||||
<LegendFont FontName="宋体" FontSize="12"/>
|
||||
|
||||
@ -0,0 +1,7 @@
|
||||
package org.meteoinfo.math.integrate;
|
||||
|
||||
public enum IntegrateMethod {
|
||||
RK45,
|
||||
RK23,
|
||||
DOP853;
|
||||
}
|
||||
@ -3,14 +3,18 @@ package org.meteoinfo.math.integrate;
|
||||
import org.apache.commons.math4.legacy.ode.ContinuousOutputModel;
|
||||
import org.apache.commons.math4.legacy.ode.FirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math4.legacy.ode.FirstOrderIntegrator;
|
||||
import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince54Integrator;
|
||||
import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince853Integrator;
|
||||
import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
|
||||
import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
|
||||
import org.meteoinfo.ndarray.Array;
|
||||
import org.meteoinfo.ndarray.DataType;
|
||||
import org.meteoinfo.ndarray.Index2D;
|
||||
import org.meteoinfo.ndarray.IndexIterator;
|
||||
|
||||
import java.io.*;
|
||||
import java.util.ArrayList;
|
||||
import java.util.List;
|
||||
|
||||
public class IntegrateUtil {
|
||||
|
||||
@ -57,4 +61,83 @@ public class IntegrateUtil {
|
||||
|
||||
return r;
|
||||
}
|
||||
|
||||
/**
|
||||
* Integrate a system of ordinary differential equations
|
||||
*
|
||||
* @param equations Computes the derivative of y at t
|
||||
* @param y0 Initial condition on y
|
||||
* @param t0 Integrate start time
|
||||
* @param t Integrate end time
|
||||
* @param tEval Times at which to store the computed solution
|
||||
* @param minStep Minimal step
|
||||
* @param maxStep Maximal step
|
||||
* @param rtol Allowed relative error
|
||||
* @param atol Allowed absolute error
|
||||
* @return Array containing the value of y for each desired time in t, with the initial value y0 in the first row.
|
||||
*/
|
||||
public static Array[] ode45(FirstOrderDifferentialEquations equations, Array y0, double t0, double t,
|
||||
Array tEval, Double minStep, Double maxStep, Double rtol, Double atol) throws IOException, ClassNotFoundException {
|
||||
y0 = y0.copyIfView();
|
||||
|
||||
int ny0 = (int) y0.getSize();
|
||||
minStep = minStep == null ? 1.0e-6 : minStep;
|
||||
maxStep = maxStep == null ? 100.0 : maxStep;
|
||||
rtol = rtol == null ? 1.0e-6 : rtol;
|
||||
atol = atol == null ? 1.0e-6 : atol;
|
||||
FirstOrderIntegrator integrator = new DormandPrince54Integrator(minStep, maxStep, atol, rtol);
|
||||
double[] y0v = (double[]) y0.getStorage();
|
||||
double[] yDot = new double[ny0];
|
||||
List<Double> tlist = new ArrayList<>();
|
||||
List<double[]> ylist = new ArrayList<>();
|
||||
|
||||
if (tEval == null) {
|
||||
integrator.addStepHandler(new StepHandler() {
|
||||
@Override
|
||||
public void init(double t0, double[] y0, double t) {
|
||||
tlist.add(t0);
|
||||
ylist.add(y0);
|
||||
//System.out.printf("Start: t0 = %.4f, y0 = %.6f%n", t0, y0[0]);
|
||||
}
|
||||
|
||||
@Override
|
||||
public void handleStep(StepInterpolator interpolator, boolean isLast) {
|
||||
double t = interpolator.getCurrentTime();
|
||||
tlist.add(t);
|
||||
double[] y = interpolator.getInterpolatedState();
|
||||
ylist.add(y.clone());
|
||||
//System.out.printf("t = %.4f, y = %.6f%n", t, y[0]);
|
||||
}
|
||||
});
|
||||
|
||||
integrator.integrate(equations, t0, y0v, t, yDot);
|
||||
} else {
|
||||
ContinuousOutputModel solution = new ContinuousOutputModel();
|
||||
integrator.addStepHandler(solution);
|
||||
|
||||
integrator.integrate(equations, t0, y0v, t, yDot);
|
||||
|
||||
for (int i = 0; i < tEval.getSize(); i++) {
|
||||
tlist.add(tEval.getDouble(i));
|
||||
solution.setInterpolatedTime(tEval.getDouble(i));
|
||||
ylist.add(solution.getInterpolatedState());
|
||||
}
|
||||
}
|
||||
|
||||
int nt = tlist.size();
|
||||
Array rt = Array.factory(DataType.DOUBLE, new int[]{nt});
|
||||
Array r = Array.factory(DataType.DOUBLE, new int[]{ny0, nt});
|
||||
Index2D index = (Index2D) r.getIndex();
|
||||
for (int i = 0; i < nt; i++) {
|
||||
rt.setDouble(i, tlist.get(i));
|
||||
double[] interpolatedY = ylist.get(i);
|
||||
index.set1(i);
|
||||
for (int j = 0; j < ny0; j++) {
|
||||
index.set0(j);
|
||||
r.setDouble(index, interpolatedY[j]);
|
||||
}
|
||||
}
|
||||
|
||||
return new Array[]{rt, r};
|
||||
}
|
||||
}
|
||||
|
||||
@ -0,0 +1,219 @@
|
||||
|
||||
package org.meteoinfo.math.integrate;
|
||||
|
||||
import org.apache.commons.math4.legacy.ode.ContinuousOutputModel;
|
||||
import org.apache.commons.math4.legacy.ode.FirstOrderDifferentialEquations;
|
||||
import org.apache.commons.math4.legacy.ode.FirstOrderIntegrator;
|
||||
import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince54Integrator;
|
||||
import org.apache.commons.math4.legacy.ode.nonstiff.DormandPrince853Integrator;
|
||||
import org.apache.commons.math4.legacy.ode.sampling.StepHandler;
|
||||
import org.apache.commons.math4.legacy.ode.sampling.StepInterpolator;
|
||||
import org.meteoinfo.ndarray.Array;
|
||||
import org.meteoinfo.ndarray.DataType;
|
||||
import org.meteoinfo.ndarray.Index2D;
|
||||
|
||||
import java.util.ArrayList;
|
||||
import java.util.List;
|
||||
import java.util.NoSuchElementException;
|
||||
|
||||
public class ODESolver {
|
||||
|
||||
private IntegrateMethod method;
|
||||
private FirstOrderIntegrator integrator;
|
||||
private FirstOrderDifferentialEquations equations;
|
||||
private Array y0;
|
||||
private double t0;
|
||||
private double tf;
|
||||
private Array tEval;
|
||||
private double minStep = 1.0e-6;
|
||||
private double maxStep = 100.0;
|
||||
private double rTol = 1.0e-6;
|
||||
private double aTol = 1.0e-6;
|
||||
private boolean denseOutput = false;
|
||||
private ContinuousOutputModel denseOutputModel;
|
||||
private Array tResult;
|
||||
private Array yResult;
|
||||
|
||||
/**
|
||||
* Constructor
|
||||
* @param method Integration method
|
||||
* @param equations Differential equations
|
||||
* @param t0 Integration start time
|
||||
* @param tf Integration end time
|
||||
* @param y0 Initial state
|
||||
*/
|
||||
public ODESolver(String method, FirstOrderDifferentialEquations equations, double t0, double tf,
|
||||
Array y0) {
|
||||
this.method = IntegrateMethod.valueOf(method.toUpperCase());
|
||||
this.equations = equations;
|
||||
this.t0 = t0;
|
||||
this.tf = tf;
|
||||
this.y0 = y0.copyIfView();
|
||||
}
|
||||
|
||||
/**
|
||||
* Set times at which to store the computed solution
|
||||
* @param value The value
|
||||
*/
|
||||
public void setTEval(Array value) {
|
||||
this.tEval = value;
|
||||
}
|
||||
|
||||
/**
|
||||
* Set first step size
|
||||
* @param value First step size
|
||||
*/
|
||||
public void setMinStep(double value) {
|
||||
this.minStep = value;
|
||||
}
|
||||
|
||||
/**
|
||||
* Set maximum allowed step size
|
||||
* @param value Maximum allowed step size
|
||||
*/
|
||||
public void setMaxStep(double value) {
|
||||
this.maxStep = value;
|
||||
}
|
||||
|
||||
/**
|
||||
* Set relative tolerance
|
||||
* @param value Relative tolerance
|
||||
*/
|
||||
public void setRTol(double value) {
|
||||
this.rTol = value;
|
||||
}
|
||||
|
||||
/**
|
||||
* Set absolute tolerance
|
||||
* @param value Absolute tolerance
|
||||
*/
|
||||
public void setATol(double value) {
|
||||
this.aTol = value;
|
||||
}
|
||||
|
||||
/**
|
||||
* Set dense output or not
|
||||
* @param value Dense output or not
|
||||
*/
|
||||
public void setDenseOutput(boolean value) {
|
||||
this.denseOutput = value;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get time result
|
||||
* @return Time result
|
||||
*/
|
||||
public Array getTResult() {
|
||||
return this.tResult;
|
||||
}
|
||||
|
||||
/**
|
||||
* Get y result
|
||||
* @return Y result
|
||||
*/
|
||||
public Array getYResult() {
|
||||
return this.yResult;
|
||||
}
|
||||
|
||||
/**
|
||||
* Solve function
|
||||
* @throws NoSuchMethodException
|
||||
*/
|
||||
public void solve() throws NoSuchMethodException {
|
||||
switch (this.method) {
|
||||
case RK45:
|
||||
this.integrator = new DormandPrince54Integrator(minStep, maxStep, aTol, rTol);
|
||||
break;
|
||||
case DOP853:
|
||||
this.integrator = new DormandPrince853Integrator(minStep, maxStep, aTol, rTol);
|
||||
break;
|
||||
default:
|
||||
throw new NoSuchMethodException("No such");
|
||||
}
|
||||
|
||||
int ny0 = (int) y0.getSize();
|
||||
double[] y0v = (double[]) y0.getStorage();
|
||||
double[] yDot = new double[ny0];
|
||||
List<Double> tlist = new ArrayList<>();
|
||||
List<double[]> ylist = new ArrayList<>();
|
||||
|
||||
if (this.denseOutput) {
|
||||
denseOutputModel = new ContinuousOutputModel();
|
||||
integrator.addStepHandler(denseOutputModel);
|
||||
|
||||
integrator.integrate(equations, t0, y0v, tf, yDot);
|
||||
} else {
|
||||
if (tEval == null) {
|
||||
integrator.addStepHandler(new StepHandler() {
|
||||
@Override
|
||||
public void init(double t0, double[] y0, double t) {
|
||||
tlist.add(t0);
|
||||
ylist.add(y0);
|
||||
//System.out.printf("Start: t0 = %.4f, y0 = %.6f%n", t0, y0[0]);
|
||||
}
|
||||
|
||||
@Override
|
||||
public void handleStep(StepInterpolator interpolator, boolean isLast) {
|
||||
double t = interpolator.getCurrentTime();
|
||||
tlist.add(t);
|
||||
double[] y = interpolator.getInterpolatedState();
|
||||
ylist.add(y.clone());
|
||||
//System.out.printf("t = %.4f, y = %.6f%n", t, y[0]);
|
||||
}
|
||||
});
|
||||
|
||||
integrator.integrate(equations, t0, y0v, tf, yDot);
|
||||
} else {
|
||||
ContinuousOutputModel denseOutputModel = new ContinuousOutputModel();
|
||||
integrator.addStepHandler(denseOutputModel);
|
||||
|
||||
integrator.integrate(equations, t0, y0v, tf, yDot);
|
||||
|
||||
for (int i = 0; i < tEval.getSize(); i++) {
|
||||
tlist.add(tEval.getDouble(i));
|
||||
denseOutputModel.setInterpolatedTime(tEval.getDouble(i));
|
||||
ylist.add(denseOutputModel.getInterpolatedState());
|
||||
}
|
||||
}
|
||||
|
||||
int nt = tlist.size();
|
||||
tResult = Array.factory(DataType.DOUBLE, new int[]{nt});
|
||||
yResult = Array.factory(DataType.DOUBLE, new int[]{ny0, nt});
|
||||
Index2D index = (Index2D) yResult.getIndex();
|
||||
for (int i = 0; i < nt; i++) {
|
||||
tResult.setDouble(i, tlist.get(i));
|
||||
double[] interpolatedY = ylist.get(i);
|
||||
index.set1(i);
|
||||
for (int j = 0; j < ny0; j++) {
|
||||
index.set0(j);
|
||||
yResult.setDouble(index, interpolatedY[j]);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/**
|
||||
* Solve by time steps - only valid with dense output
|
||||
* @param t Time steps
|
||||
* @return Values of the solution at time steps
|
||||
*/
|
||||
public Array solve(Array t) {
|
||||
t = t.copyIfView();
|
||||
int ny0 = (int) y0.getSize();
|
||||
int nt = (int) t.getSize();
|
||||
Array y = Array.factory(DataType.DOUBLE, new int[]{ny0, nt});
|
||||
Index2D index = (Index2D) y.getIndex();
|
||||
for (int i = 0; i < t.getSize(); i++) {
|
||||
denseOutputModel.setInterpolatedTime(t.getDouble(i));
|
||||
double[] v = denseOutputModel.getInterpolatedState();
|
||||
index.set1(i);
|
||||
for (int j = 0; j < ny0; j++) {
|
||||
index.set0(j);
|
||||
y.setDouble(index, v[j]);
|
||||
}
|
||||
}
|
||||
|
||||
return y;
|
||||
}
|
||||
|
||||
}
|
||||
File diff suppressed because it is too large
Load Diff
@ -0,0 +1,130 @@
|
||||
package org.meteoinfo.math.integrate.lsoda;
|
||||
|
||||
public class SupportingFunctions {
|
||||
/**
|
||||
* @param n number of elements
|
||||
* @param v value of the vector
|
||||
* @param w weights
|
||||
* @return vector norm
|
||||
*/
|
||||
public static double vmnorm(int n, double[] v, double[] w){
|
||||
double vm = 0.0;
|
||||
for (int i=1; i<=n; i++)
|
||||
vm = Math.max(vm, Math.abs(v[i])*w[i]);
|
||||
return vm;
|
||||
}
|
||||
|
||||
/**
|
||||
* full matrix norm
|
||||
* Computes the norm of a full n by n matrix,
|
||||
* stored in the array a, that is consistent with the weighted max-norm
|
||||
* on vectors, with weights stored in the array w.
|
||||
* fnorm = max(i=1,...,n) ( w[i] * sum(j=1,...,n) fabs( a[i][j] ) / w[j] )
|
||||
* @param n dimensions
|
||||
* @param a the matrix
|
||||
* @param w weights
|
||||
* @return matrix norm
|
||||
*/
|
||||
public static double fnorm(int n, double[][] a, double[] w){
|
||||
double an = 0.0, sum;
|
||||
for (int i=1; i <= n ; i++){
|
||||
sum = 0;
|
||||
for(int j=1; j<=n; j++)
|
||||
sum += Math.abs(a[i][j])/w[j];
|
||||
an = Math.max(an, sum*w[i]);
|
||||
}
|
||||
return an;
|
||||
}
|
||||
|
||||
/**
|
||||
* compute EWT
|
||||
* @param ycur current value of y
|
||||
*/
|
||||
public static double[] ewset(double[] ycur, int itol, double[] rtol, double[] atol, int n){
|
||||
int i;
|
||||
double[] EWT = new double[n+1];
|
||||
switch (itol){
|
||||
case 1:
|
||||
for (i=1; i<=n; i++) // both are scalars
|
||||
EWT[i] = rtol[1] * Math.abs(ycur[i])+atol[1];
|
||||
break;
|
||||
case 2:
|
||||
for (i=1;i<=n; i++) //scalar rtol and array atol
|
||||
EWT[i] = rtol[1] * Math.abs(ycur[i]) + atol[i];
|
||||
break;
|
||||
case 3:
|
||||
for (i=1;i<=n; i++) //array rtol and scalar atol
|
||||
EWT[i] = rtol[i] * Math.abs(ycur[i]) + atol[1];
|
||||
break;
|
||||
case 4:
|
||||
for (i=1;i<=n; i++) // both are arrays
|
||||
EWT[i] = rtol[i] * Math.abs(ycur[i]) + atol[i];
|
||||
break;
|
||||
}
|
||||
return EWT;
|
||||
}
|
||||
|
||||
/**
|
||||
* Manages the solution of the linear system arising from
|
||||
* a chord iteration. It is called if miter != 0.
|
||||
* @param y the right-hand side vector on input, and the solution vector on output.
|
||||
* @return solution of the linear system
|
||||
*/
|
||||
public static double[] solsy(double[] y, double[][] wm, int n, int[] ipvt, int miter){
|
||||
if (miter != 2){
|
||||
System.out.println("solsy: miter != 2");
|
||||
return y;
|
||||
}
|
||||
return Utility.solveLinearSys(wm,n,ipvt,y);
|
||||
}
|
||||
|
||||
/**
|
||||
* interpolate at the output station
|
||||
* @param t value of independent variable
|
||||
* @param k integer that specifies the desired derivative order
|
||||
*/
|
||||
public static double[] intdy(double t, int k, double[] dky, int nq, double tn,
|
||||
double hu, double ETA, double h, double[][] yh, int n){
|
||||
int ic ,jp1;
|
||||
double r, s, tp;
|
||||
|
||||
if (k < 0 || k > nq){
|
||||
System.err.printf("intdy: k = %d illegal", k);
|
||||
dky[0] = -1;
|
||||
return dky;
|
||||
}
|
||||
tp = tn - hu -100.0 * ETA * (tn + hu);
|
||||
if ((t-tp) * (t-tn)> 0.0){
|
||||
System.err.printf("intdy: t = %f illegal\n" +
|
||||
"t not in interval tcur-hu to tcur\n", t);
|
||||
dky[0] = -2;
|
||||
return dky;
|
||||
}
|
||||
|
||||
s = (t-tn)/h;
|
||||
ic = 1;
|
||||
for (int i = nq+1-k; i <= nq; i++){
|
||||
ic *= i; // ic = nq(nq-1)...(nq-k+1)
|
||||
}
|
||||
for (int i=1; i<=n; i++)
|
||||
dky[i] = ic* yh[nq+1][i];
|
||||
|
||||
// use Horner's rule
|
||||
for (int j = nq-1; j >= k; j--){
|
||||
jp1 = j+1; // order of current derivatives
|
||||
ic = 1;
|
||||
for (int jj = jp1-k; jj<=j;jj++)
|
||||
ic *= jj; // coefficient
|
||||
for (int i=1; i<=n; i++)
|
||||
dky[i] = ic * yh[jp1][i] + s * dky[i];
|
||||
}
|
||||
if (k == 0){
|
||||
return dky;
|
||||
}
|
||||
r = Math.pow(h, -k);
|
||||
for (int i=1; i<=n; i++)
|
||||
dky[i] = r*dky[i];
|
||||
|
||||
return dky;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,314 @@
|
||||
package org.meteoinfo.math.integrate.lsoda;
|
||||
|
||||
import java.util.Arrays;
|
||||
|
||||
public class Utility {
|
||||
/**
|
||||
* ddot.c
|
||||
* @param n number of elements
|
||||
* @param dx n+1 vector (dx[0] is not used)
|
||||
* @param incx storage spacing between elements of dx
|
||||
* @param dy n+1 vector (dy[0] is not used)
|
||||
* @param incy storage spacing between elements of dy
|
||||
* @return dot product of dx and dy
|
||||
*/
|
||||
public static double Dot(int n, double[] dx, int incx, double[] dy, int incy){
|
||||
double dotProd = 0; // result
|
||||
int ix = 1, iy = 1; // index of dx and dy
|
||||
if (n <= 0)
|
||||
return dotProd;
|
||||
|
||||
// unequal or non-positive increments
|
||||
if( incx != incy || incx < 1){
|
||||
|
||||
if (incx < 0)
|
||||
ix = (-n+1) * incx + 1;
|
||||
if (incx < 0)
|
||||
iy = (-n+1) * incy + 1;
|
||||
|
||||
for(int i=1;i<=n;i++){
|
||||
dotProd += dx[ix] * dy[iy];
|
||||
ix += incx;
|
||||
iy += incy;
|
||||
}
|
||||
return dotProd;
|
||||
}
|
||||
|
||||
// Both increments equal to 1
|
||||
if(incx == 1){
|
||||
int m = n % 5;
|
||||
// Clean-up loop so remaining vector length is a multiple of 5
|
||||
if (m != 0){
|
||||
for (int i=1; i<=m;i++)
|
||||
dotProd += dx[i] * dy[i];
|
||||
if (n < 5)
|
||||
return dotProd;
|
||||
}
|
||||
for (int i=m+1; i<=n;i=i+5)
|
||||
dotProd += dx[i] * dy[i] + dx[i+1] * dy[i+1] +
|
||||
dx[i+2] * dy[i+2] + dx[i+3] * dy[i+3] +
|
||||
dx[i+4] * dy[i+4];
|
||||
return dotProd;
|
||||
}
|
||||
// positive equal non-unit increments
|
||||
for (int i=1; i<= n*incx; i=i+incx)
|
||||
dotProd += dx[i] * dy[i];
|
||||
|
||||
return dotProd;
|
||||
}
|
||||
|
||||
/**
|
||||
* idamax.c
|
||||
* Find largest component of double vector dx
|
||||
* @param n number of elements
|
||||
* @param dx n+1 vector (dx[0] is not used)
|
||||
* @param incx storage spacing between elements of dx
|
||||
* @return smallest index
|
||||
*/
|
||||
public static int findMaxMagnitude(int n, double[] dx, int incx){
|
||||
double maxValue, curValue;
|
||||
int xindex = 0;
|
||||
|
||||
if (n <= 0)
|
||||
return xindex;
|
||||
xindex = 1;
|
||||
if (n <= 1 || incx <= 0)
|
||||
return xindex;
|
||||
|
||||
maxValue = Math.abs(dx[1]);
|
||||
// increments are not euqal to 1
|
||||
if (incx != 1){
|
||||
int curIndex = 2;
|
||||
for (int i=1+incx; i <= n*incx; i=i+incx){
|
||||
curValue = Math.abs(dx[i]);
|
||||
if (curValue > maxValue){
|
||||
xindex = curIndex;
|
||||
maxValue = curValue;
|
||||
}
|
||||
curIndex++;
|
||||
}
|
||||
return xindex;
|
||||
}
|
||||
|
||||
// increments are equal to 1.
|
||||
for (int i=2; i <= n; i++){
|
||||
curValue = Math.abs(dx[i]);
|
||||
if (curValue > maxValue){
|
||||
xindex = i;
|
||||
maxValue = curValue;
|
||||
}
|
||||
}
|
||||
return xindex;
|
||||
}
|
||||
|
||||
/**
|
||||
* dscal.c
|
||||
* scalar vector multiplication
|
||||
* @param n number of elements
|
||||
* @param da double scalar multiplier
|
||||
* @param dx start+n-1 vector (dx[0] is not used)
|
||||
* @param incx storage spacing between elements of dx (should be positive)
|
||||
* @param start start position
|
||||
* @return da*dx
|
||||
*/
|
||||
public static double[] calAX(int n, double da, double[] dx, int incx, int start){
|
||||
if (n <= 0 || incx * (n-1) + start >= dx.length){
|
||||
System.out.println("function calAX: illegal input");
|
||||
return dx;
|
||||
}
|
||||
|
||||
// increments are not equal to 1
|
||||
if (incx != 1){
|
||||
for (int i = start; i <= (n-1)*incx + start; i = i+incx)
|
||||
dx[i] = da * dx[i];
|
||||
return dx;
|
||||
}
|
||||
|
||||
// increments are equal to 1
|
||||
int m = n % 5;
|
||||
if ( m != 0 ) {
|
||||
for (int i = start ; i <= m+start - 1 ; i++ )
|
||||
dx[i] = da * dx[i];
|
||||
if ( n < 5 )
|
||||
return dx;
|
||||
}
|
||||
for (int i = m + start; i <= n + start - 1; i = i + 5 ) {
|
||||
dx[i] = da * dx[i];
|
||||
dx[i+1] = da * dx[i+1];
|
||||
dx[i+2] = da * dx[i+2];
|
||||
dx[i+3] = da * dx[i+3];
|
||||
dx[i+4] = da * dx[i+4];
|
||||
}
|
||||
return dx;
|
||||
}
|
||||
|
||||
/**
|
||||
* daxpy.c
|
||||
* @param n number of elements
|
||||
* @param da double scalar multiplier
|
||||
* @param dx start+n-1 vector (dx[0] is not used)
|
||||
* @param incx storage spacing between elements of dx
|
||||
* @param dy start+n-1 vector (dy[0] is not used)
|
||||
* @param incy storage spacing between elements of dy
|
||||
* @param ystart start position of y
|
||||
* @param xstart start position of x
|
||||
* @return da * dx + dy
|
||||
*/
|
||||
public static double[] calAXPlusY(int n, double da, double[] dx, int incx, double[] dy, int incy, int xstart, int ystart){
|
||||
int ix, iy; // index
|
||||
|
||||
if (incx * (n-1) + xstart >= dx.length||incy * (n-1) + ystart >= dy.length||n<=0){
|
||||
System.out.println("function calAXPlusY: illegal input");
|
||||
return dy;
|
||||
}
|
||||
|
||||
if ( da == 0.0)
|
||||
return dy;
|
||||
|
||||
// nonequal or nonpositive increments
|
||||
if (incx != incy || incx < 1){
|
||||
ix = xstart;
|
||||
iy = ystart;
|
||||
if (incx < 0)
|
||||
ix = (-n+1) * incx + xstart;
|
||||
if (incx < 0)
|
||||
iy = (-n+1) * incy + ystart;
|
||||
for (int i=1; i<=n; i++){
|
||||
dy[iy] += da * dx[ix];
|
||||
ix += incx;
|
||||
iy += incy;
|
||||
}
|
||||
return dy;
|
||||
}
|
||||
|
||||
// increments are equal to 1
|
||||
if (incx == 1){
|
||||
// Clean-up loop
|
||||
int m = n % 4;
|
||||
if (m != 0){
|
||||
for (int i=1; i<=m; i++)
|
||||
dy[i+ystart-1] += da * dx[i+xstart-1];
|
||||
if (n < 4)
|
||||
return dy;
|
||||
}
|
||||
for (int i=m+1; i <= n; i=i+4){
|
||||
dy[i+ystart-1] = dy[i+ystart-1] + da * dx[i+xstart-1];
|
||||
dy[i+ystart] = dy[i+ystart] + da * dx[i+xstart];
|
||||
dy[i+ystart+1] = dy[i+ystart+1] + da * dx[i+xstart+1];
|
||||
dy[i+ystart+2] = dy[i+ystart+2] + da * dx[i+xstart+2];
|
||||
}
|
||||
return dy;
|
||||
}
|
||||
|
||||
// incrememts are equal, positive but not unit
|
||||
for (int i = 1; i <= n*incx; i=i+incx)
|
||||
dy[ystart+i-1] += da*dx[xstart+i-1];
|
||||
return dy;
|
||||
}
|
||||
|
||||
|
||||
/**
|
||||
* LU decomposition using partial pivoting
|
||||
* @param a (n+1)x(n+1) matrix
|
||||
* @param n dimension
|
||||
* @return a, ipvt, info
|
||||
* if info != 0, the matrix is singular
|
||||
*/
|
||||
public static ReturningValues LUDecomposition(double[][] a, int n){
|
||||
double t;
|
||||
int info = 0;
|
||||
int[] ipvt = new int[n+1];
|
||||
|
||||
for (int k=1 ; k<=n-1 ; k++){
|
||||
// Find j = pivot index
|
||||
int j = findMaxMagnitude(n-k+1, Arrays.copyOfRange(a[k],k-1,n+1),1)+k-1;
|
||||
ipvt[k] = j;
|
||||
|
||||
// zero pivot implies this row already triangularised
|
||||
if (a[k][j] == 0.0){
|
||||
info = k;
|
||||
continue;
|
||||
}
|
||||
|
||||
// interchange
|
||||
if(j != k){
|
||||
t = a[k][j];
|
||||
a[k][j] = a[k][k];
|
||||
a[k][k] = t;
|
||||
}
|
||||
|
||||
// compute multipliers
|
||||
t = - 1.0/ a[k][k];
|
||||
a[k] = calAX(n-k, t, a[k],1, k+1);
|
||||
|
||||
// column elimination with row indexing
|
||||
for (int i=k+1; i<=n;i++){
|
||||
t = a[i][j];
|
||||
if(j != k){
|
||||
a[i][j] = a[i][k];
|
||||
a[i][k] = t;
|
||||
}
|
||||
a[i] = calAXPlusY(n-k, t, a[k], 1, a[i], 1,k+1,k+1);
|
||||
}
|
||||
}
|
||||
|
||||
ipvt[n] = n;
|
||||
if(a[n][n] == 0.0)
|
||||
info = n;
|
||||
|
||||
return new ReturningValues(a, ipvt, info);
|
||||
}
|
||||
|
||||
/**
|
||||
* solves the linear system using the factors computed by LUDecomposition
|
||||
* @param a (n+1)x(n+1) matrix
|
||||
* @param n row dimension of a
|
||||
* @param ipvt the pivot vector from LUDecomposition
|
||||
* @param b the right hand side vector
|
||||
* @return b the solution vector x.
|
||||
*/
|
||||
public static double[] solveLinearSys(double[][] a, int n, int[] ipvt, double[] b){
|
||||
int j;
|
||||
double t;
|
||||
|
||||
// Ax = LUx = b; Let Ux = y, so Ly = b.
|
||||
// 1. Solve L * y = b.
|
||||
for(int k=1; k <= n; k++){
|
||||
t = Dot(k-1, Arrays.copyOfRange(a[k],0,k),1,b,1);
|
||||
b[k] = (b[k]-t) / a[k][k];
|
||||
}
|
||||
|
||||
// 2. Solve U * x = y.
|
||||
for (int k=n-1; k >= 1; k--){
|
||||
b[k] = b[k] + Dot(n-k, Arrays.copyOfRange(a[k],k,n+1),1, Arrays.copyOfRange(b,k,n+1),1);
|
||||
j = ipvt[k];
|
||||
if (j != k){
|
||||
t = b[j];
|
||||
b[j] = b[k];
|
||||
b[k] = t;
|
||||
}
|
||||
}
|
||||
return b;
|
||||
}
|
||||
|
||||
public static Double[] transformYVec(double[] y){
|
||||
int len = y.length-1;
|
||||
Double[] yvec = new Double[len];
|
||||
for (int i=0; i<len;i++)
|
||||
yvec[i] = y[i+1];
|
||||
return yvec;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
|
||||
class ReturningValues{
|
||||
public final double[][] a;
|
||||
public final int info;
|
||||
public final int[] ipvt;
|
||||
public ReturningValues(double[][] a, int[] ipvt, int info){
|
||||
this.info = info;
|
||||
this.a = a;
|
||||
this.ipvt = ipvt;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,15 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class EwtException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public EwtException(Number i, Number value){
|
||||
super();
|
||||
errorMsg ="lsoda: ewt["+i+"] = "+value+" <= 0.0";
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getMessage() {
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,15 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class ExceedMaxStepsException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public ExceedMaxStepsException(int mxstep){
|
||||
super();
|
||||
errorMsg = "lsoda: "+mxstep+" steps taken before reaching tout";
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getMessage() {
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,14 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class IllegalInputException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public IllegalInputException(String param, Number value){
|
||||
super();
|
||||
errorMsg = "lsoda: illegal "+param+"="+value;
|
||||
}
|
||||
@Override
|
||||
public String getMessage() {
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,30 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class IllegalTException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public IllegalTException(String type, double t, double tout, Number tt){
|
||||
switch (type){
|
||||
case "tout_behind_t":
|
||||
errorMsg = "lsoda: tout ="+tout+" behind t = "+t +
|
||||
"\n integration direction is given by "+tt;
|
||||
break;
|
||||
case "tcrit_behind_tout":
|
||||
errorMsg = "lsoda: itask = 4 or 5 and tcrit behind tout";
|
||||
break;
|
||||
case "tcrit_behind_tcur":
|
||||
errorMsg = "lsoda: itask = 4 or 5 and tcrit behind tcur";
|
||||
break;
|
||||
case "tout_close_to_t":
|
||||
errorMsg = "lsoda: tout too close to t to start integration";
|
||||
break;
|
||||
case "tout_behind_tcur_hu":
|
||||
errorMsg = "lsoda: itask = "+tt+" and tout behind tcur - hu";
|
||||
break;
|
||||
}
|
||||
}
|
||||
@Override
|
||||
public String getMessage(){
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,16 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class InterpolationException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public InterpolationException(int itask, double tout){
|
||||
super();
|
||||
errorMsg = "lsoda: trouble from intdy, itask = "+itask+
|
||||
", tout = "+tout;
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getMessage() {
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,21 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class IstateException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public IstateException(int istate){
|
||||
super();
|
||||
errorMsg = "lsoda: illegal istate ="+istate;
|
||||
}
|
||||
|
||||
public IstateException(int istate, int init){
|
||||
super();
|
||||
if (istate>1 && init==0)
|
||||
errorMsg = "lsoda: istate > 1 but lsoda not initialised";
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getMessage() {
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,20 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class NeqException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public NeqException(int neq){
|
||||
super();
|
||||
errorMsg = "lsoda: neq = "+neq+" is less than 1";
|
||||
}
|
||||
|
||||
public NeqException(int neq, int istate){
|
||||
super();
|
||||
errorMsg = "lsoda: istate = "+istate+" and neq increased";
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getMessage() {
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,16 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class RepeatedInputException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public RepeatedInputException(){
|
||||
super();
|
||||
errorMsg ="lsoda: repeated calls with istate = 1 and tout = t\n"+
|
||||
" run aborted.. apparent infinite loop";
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getMessage() {
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,23 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class TestsFailException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public TestsFailException(String type, double tn, double h){
|
||||
super();
|
||||
errorMsg = "lsoda: at t =" + tn + " and step size h = " + h+"\n";
|
||||
switch (type){
|
||||
case "error_test":
|
||||
errorMsg += " error test failed repeatedly or with abs(h)=hmin";
|
||||
break;
|
||||
case "convergence_test":
|
||||
errorMsg += " corrector convergence failed repeatedly or with abs(h)=hmin";
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getMessage() {
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,23 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.exception;
|
||||
|
||||
public class TooMuchAccuracyException extends RuntimeException{
|
||||
private String errorMsg;
|
||||
|
||||
public TooMuchAccuracyException(double tolsf){
|
||||
super();
|
||||
errorMsg = "lsoda: at start of problem, too much accuracy\n" +
|
||||
" requested for precision of machine,\n" +
|
||||
" suggested scaling factor = "+tolsf;
|
||||
}
|
||||
|
||||
public TooMuchAccuracyException(double t, double tolsf){
|
||||
errorMsg = "lsoda: at t = " + t + ", too much accuracy requested\n" +
|
||||
" for precision of machine, suggested\n" +
|
||||
" scaling factor = " + tolsf;
|
||||
}
|
||||
|
||||
@Override
|
||||
public String getMessage() {
|
||||
return errorMsg;
|
||||
}
|
||||
}
|
||||
@ -0,0 +1,42 @@
|
||||
package org.meteoinfo.math.integrate.lsoda.tools;
|
||||
|
||||
|
||||
import java.io.BufferedWriter;
|
||||
import java.io.File;
|
||||
import java.io.FileWriter;
|
||||
import java.io.IOException;
|
||||
|
||||
public class Data {
|
||||
private String filePath;
|
||||
private BufferedWriter bufferedWriter;
|
||||
private int n;
|
||||
|
||||
public Data(String filepath, String fileName, int neq) throws IOException {
|
||||
n=neq;
|
||||
filePath = filepath + "/" + fileName;
|
||||
// header
|
||||
bufferedWriter = new BufferedWriter(new FileWriter(filePath));
|
||||
bufferedWriter.write("t,");
|
||||
for (int i=1; i<neq; i++)
|
||||
bufferedWriter.write("y_"+i+",");
|
||||
bufferedWriter.write("y_"+neq);
|
||||
bufferedWriter.flush();
|
||||
}
|
||||
|
||||
public void write(double t, double[] y) {
|
||||
try{
|
||||
bufferedWriter.newLine();
|
||||
bufferedWriter.write(t +",");
|
||||
for (int i=1; i<n; i++)
|
||||
bufferedWriter.write(y[i]+",");
|
||||
bufferedWriter.write(Double.toString(y[n]));
|
||||
bufferedWriter.flush();
|
||||
} catch (IOException e) {
|
||||
System.out.println("Error occurred in file writing.");
|
||||
}
|
||||
}
|
||||
|
||||
public void closeWriter() throws IOException {
|
||||
bufferedWriter.close();
|
||||
}
|
||||
}
|
||||
3
pom.xml
3
pom.xml
@ -30,12 +30,13 @@
|
||||
<module>meteoinfo-projection</module>
|
||||
<module>meteoinfo-mkl</module>
|
||||
<module>meteoinfo-render2d</module>
|
||||
<module>meteoinfo-jython</module>
|
||||
</modules>
|
||||
|
||||
<properties>
|
||||
<project.build.sourceEncoding>UTF-8</project.build.sourceEncoding>
|
||||
<java.version>1.8</java.version>
|
||||
<revision>4.0.4</revision>
|
||||
<revision>4.0.5</revision>
|
||||
<maven.compiler.source>8</maven.compiler.source>
|
||||
<maven.compiler.target>8</maven.compiler.target>
|
||||
<maven.compiler.release>8</maven.compiler.release>
|
||||
|
||||
Loading…
x
Reference in New Issue
Block a user