acb_ode, numerical solution of linear ODEs with polynomial coefficients

Makefile.in

@@ -217,7 +217,7 @@ HEADER_DIRS	:=	\
 		double_interval	dlog								\
 		arb_fmpz_poly	arb_fpwrap							\
 		acb_dft			acb_elliptic	acb_modular			\
-		acb_dirichlet	acb_theta							\
+		acb_dirichlet	acb_theta		acb_ode				\
 		dirichlet		bernoulli		hypgeom				\
 		bool_mat		partitions							\
 															\

dev/check_examples.sh

@@ -169,6 +169,9 @@ then
 elif test "$1" = "hilbert_matrix_ca";
 then
     echo "hilbert_matrix_ca....SKIPPED"
+elif test "$1" = "ode";
+then
+    echo "ode....SKIPPED"
 elif test "$1" = "huge_expr";
 then
     echo "huge_expr....SKIPPED"

doc/source/acb.rst

@@ -1206,6 +1206,14 @@ Vector functions
 
     Swaps the entries of *vec1* and *vec2*.
 
+.. function:: void _acb_vec_get_mag(mag_t bound, acb_srcptr vec, slong len)
+
+    Sets *bound* to an upper bound for the entries in *vec*.
+
+.. function:: void _acb_vec_get_mag_lower(mag_t bound, acb_srcptr vec, slong len)
+
+    Sets *bound* to an lower bound for the absolute values of the entries in *vec*.
+
 .. function:: void _acb_vec_get_real(arb_ptr re, acb_srcptr vec, slong len)
 
 .. function:: void _acb_vec_get_imag(arb_ptr im, acb_srcptr vec, slong len)

doc/source/acb_poly.rst

@@ -1205,3 +1205,25 @@ Root-finding
     If this function returns zero, then the signs of the imaginary parts
     are not known for certain, based on the accuracy of the inputs
     and the working precision *prec*.
+
+
+Vector functions
+-------------------------------------------------------------------------------
+
+.. function:: acb_poly_struct * _acb_poly_vec_init(slong n)
+
+.. function:: void _acb_poly_vec_clear(acb_poly_struct * vec, slong n)
+
+.. function:: void _acb_poly_vec_set(acb_poly_struct * dest, const acb_poly_struct * src, slong n)
+
+.. function:: void _acb_poly_vec_set_block(acb_poly_struct * dest, const acb_poly_struct * src, slong n, slong base, slong len)
+
+.. function:: slong _acb_poly_vec_length(const acb_poly_struct * vec, slong n)
+
+.. function:: void _acb_poly_vec_fit_length(acb_poly_struct * vec, slong n, slong len)
+
+.. function:: void _acb_poly_vec_set_length(acb_poly_struct * vec, slong n, slong len)
+
+.. function:: void _acb_poly_vec_normalise(acb_poly_struct * vec, slong n)
+
+.. function:: int _acb_poly_vec_overlaps(acb_poly_struct * vec1, acb_poly_struct * vec2, slong n)

doc/source/index.rst

@@ -209,6 +209,7 @@ Real and complex numbers
    acb_modular.rst
    acb_theta.rst
    acb_dirichlet.rst
+   acb_ode.rst
    bernoulli.rst
    hypgeom.rst
    partitions.rst

examples/ode_basic.c

@@ -0,0 +1,69 @@
+#include <math.h>
+
+#include "acb_types.h"
+#include "acb.h"
+#include "acb_mat.h"
+#include "acb_ode.h"
+#include "gr.h"
+#include "gr_ore_poly.h"
+
+int main(void)
+{
+    gr_ctx_t ZZ, Pol, Dop;
+    gr_ptr dop;
+
+    int status = GR_SUCCESS;
+
+    slong prec = 50;
+
+    // Differential operator
+
+    // XXX fmpz_poly? fmpq_poly??
+    // XXX std derivative
+    gr_ctx_init_fmpz(ZZ);
+    gr_ctx_init_gr_poly(Pol, ZZ);
+    gr_ctx_init_gr_ore_poly(Dop, Pol, 0, ORE_ALGEBRA_EULER_DERIVATIVE);
+    status |= gr_ctx_set_gen_name(Pol, "z");
+    status |= gr_ctx_set_gen_name(Dop, "Tz");
+    GR_TMP_INIT(dop, Dop);
+    status |= gr_ore_poly_set_str(dop,  "(z^2 + 1)*Tz^2 + (z^2 - 1)*Tz", Dop);
+    // status |= gr_ore_poly_set_str(dop,  "4*Tz^2 - 4*Tz - z^2 + 8*z - 11", Dop);
+    // status |= gr_ore_poly_set_str(dop, "Tz^6 - 6*Tz^5 + 12*Tz^4 - 10*Tz^3 + 3*Tz^2 + z^2", Dop);
+    // status |= gr_ore_poly_set_str(dop, "4*Tz^2 - 4*Tz - z^2 + 8*z - 11", Dop);
+    GR_MUST_SUCCEED(status);
+
+    // flint_printf("dop = %{gr}\n", dop, Dop);
+
+    // Evaluation point
+
+    acb_t pt;
+    acb_init(pt);
+    // acb_set_d(pt, .5);
+    acb_const_pi(pt, prec + 4);
+    arb_one(acb_imagref(pt));
+    acb_inv(pt, pt, prec + 4);
+
+    // Output matrix
+
+    acb_mat_t mat;
+    slong dop_order = gr_ore_poly_length(dop, Dop) - 1;
+    acb_mat_init(mat, dop_order, dop_order);
+
+    // Numerical solution
+
+    status |= acb_ode_fundamental_matrix(mat, dop, Dop, NULL, NULL, pt, 0, prec);
+    GR_MUST_SUCCEED(status);
+
+    flint_printf("\n------\n");
+    acb_mat_printd(mat, prec+4);
+    flint_printf("\n%f\n", prec/log2(10.));
+
+    acb_clear(pt);
+    acb_mat_clear(mat);
+    GR_TMP_CLEAR(dop, Dop);
+
+    gr_ctx_clear(Dop);
+    gr_ctx_clear(Pol);
+
+    flint_cleanup_master();
+}

examples/ode_fundamental_matrix_manual_exponents.c

@@ -0,0 +1,144 @@
+#include "acb_types.h"
+#include "acb.h"
+#include "acb_mat.h"
+#include "acb_ode.h"
+#include "gr.h"
+#include "gr_ore_poly.h"
+
+// XXX also manual sing?
+
+void
+fundamental_matrix(const char * dop_str,
+                   const acb_ode_exponents_struct * expos,
+                   double pt_d)
+{
+    gr_ctx_t CC, Pol, Dop;
+    gr_ptr dop;
+
+    slong prec = 30;
+
+    int status = GR_SUCCESS;
+
+    gr_ctx_init_complex_acb(CC, prec);
+    gr_ctx_init_gr_poly(Pol, CC);
+    gr_ctx_init_gr_ore_poly(Dop, Pol, 0, ORE_ALGEBRA_EULER_DERIVATIVE);
+
+    GR_TMP_INIT(dop, Dop);
+
+    status |= gr_ctx_set_gen_name(Pol, "z");
+    status |= gr_ctx_set_gen_name(Dop, "Tz");
+    status |= gr_ore_poly_set_str(dop, dop_str, Dop);
+
+    status |= gr_println(dop, Dop);
+    flint_printf("\n");
+
+    GR_MUST_SUCCEED(status);
+
+    slong dop_order = gr_ore_poly_length(dop, Dop) - 1;
+
+    acb_mat_t mat;
+    acb_mat_init(mat, dop_order, dop_order);
+
+    acb_t pt;
+    acb_init(pt);
+    acb_set_d(pt, pt_d);
+
+    GR_MUST_SUCCEED(acb_ode_fundamental_matrix(mat, dop, Dop, expos, NULL, pt, 0, prec));
+
+    flint_printf("%{acb_mat}\n\n", mat);
+    flint_printf("--------\n\n");
+
+    acb_mat_clear(mat);
+    GR_TMP_CLEAR(dop, Dop);
+    gr_ctx_clear(Dop);
+    gr_ctx_clear(Pol);
+    gr_ctx_clear(CC);
+    acb_clear(pt);
+}
+
+
+void
+apery(void)
+{
+    acb_ode_shift_struct shift[1] = {{ .n = 0, .mult = 3 }};
+    acb_ode_group_struct grp[1] = {{ .nshifts = 1, .shifts = shift }};
+    acb_init(grp->leader);
+    acb_zero(grp->leader);
+    acb_ode_exponents_struct expos[1] = {{ .ngroups = 1, .groups = grp }};
+
+    fundamental_matrix(
+            "(z^2 - 34*z + 1)*Tz^3 + (3*z^2 - 51*z)*Tz^2 + (3*z^2 - 27*z)*Tz + z^2 - 5*z",
+            expos,
+            0.015625);
+
+    acb_clear(grp->leader);
+}
+
+
+void
+multiple_shifts(void)
+{
+    const char * dop = "Tz^6 - 6*Tz^5 + 12*Tz^4 - 10*Tz^3 + 3*Tz^2 + z^2";
+
+    acb_ode_shift_struct shift[3] = {
+        { .n = 0, .mult = 2 },
+        { .n = 1, .mult = 3 },
+        { .n = 3, .mult = 1 },
+    };
+    acb_ode_group_struct grp[1] = {{ .nshifts = 3, .shifts = shift }};
+    acb_init(grp->leader);
+    acb_zero(grp->leader);
+    acb_ode_exponents_struct expos[1] = {{ .ngroups = 1, .groups = grp }};
+
+    fundamental_matrix(dop, expos, 2.);
+
+    acb_clear(grp->leader);
+}
+
+
+void
+whittaker(void)
+{
+    slong prec = 30;
+
+    acb_t kappa, mu, half;
+    acb_init(kappa);
+    acb_init(mu);
+    acb_init(half);
+
+    const char * dop = "4*Tz^2 - 4*Tz - z^2 + 8*z - 11";
+    acb_set_si(kappa, 2);
+    acb_set_si(mu, 3);
+    acb_sqrt(mu, mu, prec);
+    acb_set_d(half, .5);
+
+    acb_ode_shift_struct shift[1] = { { .n = 0, .mult = 1 } };
+    acb_ode_group_struct grp[2] = {
+        { .nshifts = 1, .shifts = shift },
+        { .nshifts = 1, .shifts = shift },
+    };
+    acb_init(grp[0].leader);
+    acb_sub(grp[0].leader, half, mu, prec);
+    acb_init(grp[1].leader);
+    acb_add(grp[1].leader, half, mu, prec);
+    acb_ode_exponents_struct expos[1] = {{ .ngroups = 2, .groups = grp }};
+
+    fundamental_matrix(dop, expos, 2.);
+
+    acb_clear(grp[0].leader);
+    acb_clear(grp[1].leader);
+    acb_clear(kappa);
+    acb_clear(mu);
+    acb_clear(half);
+}
+
+
+int
+main(void)
+{
+    apery();
+    multiple_shifts();
+    whittaker();
+
+    flint_cleanup_master();
+}