Skip to content

Commit a8f91a7

Browse files
mezzarobbamezzarobba
committed
[WIP] Solutions of linear ODEs with coeffs in ℂ[x]
1 parent 65198c1 commit a8f91a7

49 files changed

Lines changed: 5820 additions & 1 deletion

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

‎Makefile.in‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -217,7 +217,7 @@ HEADER_DIRS := \
217217
double_interval dlog \
218218
arb_fmpz_poly arb_fpwrap \
219219
acb_dft acb_elliptic acb_modular \
220-
acb_dirichlet acb_theta \
220+
acb_dirichlet acb_theta acb_ode \
221221
dirichlet bernoulli hypgeom \
222222
bool_mat partitions \
223223
\

‎dev/check_examples.sh‎

Lines changed: 3 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -169,6 +169,9 @@ then
169169
elif test "$1" = "hilbert_matrix_ca";
170170
then
171171
echo "hilbert_matrix_ca....SKIPPED"
172+
elif test "$1" = "ode";
173+
then
174+
echo "ode....SKIPPED"
172175
elif test "$1" = "huge_expr";
173176
then
174177
echo "huge_expr....SKIPPED"

‎doc/source/acb.rst‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1206,6 +1206,14 @@ Vector functions
12061206

12071207
Swaps the entries of *vec1* and *vec2*.
12081208

1209+
.. function:: void _acb_vec_get_mag(mag_t bound, acb_srcptr vec, slong len)
1210+
1211+
Sets *bound* to an upper bound for the entries in *vec*.
1212+
1213+
.. function:: void _acb_vec_get_mag_lower(mag_t bound, acb_srcptr vec, slong len)
1214+
1215+
Sets *bound* to an lower bound for the absolute values of the entries in *vec*.
1216+
12091217
.. function:: void _acb_vec_get_real(arb_ptr re, acb_srcptr vec, slong len)
12101218

12111219
.. function:: void _acb_vec_get_imag(arb_ptr im, acb_srcptr vec, slong len)

‎doc/source/acb_ode.rst‎

Lines changed: 654 additions & 0 deletions
Large diffs are not rendered by default.

‎doc/source/acb_poly.rst‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1205,3 +1205,25 @@ Root-finding
12051205
If this function returns zero, then the signs of the imaginary parts
12061206
are not known for certain, based on the accuracy of the inputs
12071207
and the working precision *prec*.
1208+
1209+
1210+
Vector functions
1211+
-------------------------------------------------------------------------------
1212+
1213+
.. function:: acb_poly_struct * _acb_poly_vec_init(slong n)
1214+
1215+
.. function:: void _acb_poly_vec_clear(acb_poly_struct * vec, slong n)
1216+
1217+
.. function:: void _acb_poly_vec_set(acb_poly_struct * dest, const acb_poly_struct * src, slong n)
1218+
1219+
.. function:: void _acb_poly_vec_set_block(acb_poly_struct * dest, const acb_poly_struct * src, slong n, slong base, slong len)
1220+
1221+
.. function:: slong _acb_poly_vec_length(const acb_poly_struct * vec, slong n)
1222+
1223+
.. function:: void _acb_poly_vec_fit_length(acb_poly_struct * vec, slong n, slong len)
1224+
1225+
.. function:: void _acb_poly_vec_set_length(acb_poly_struct * vec, slong n, slong len)
1226+
1227+
.. function:: void _acb_poly_vec_normalise(acb_poly_struct * vec, slong n)
1228+
1229+
.. function:: int _acb_poly_vec_overlaps(acb_poly_struct * vec1, acb_poly_struct * vec2, slong n)

‎doc/source/index.rst‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -209,6 +209,7 @@ Real and complex numbers
209209
acb_modular.rst
210210
acb_theta.rst
211211
acb_dirichlet.rst
212+
acb_ode.rst
212213
bernoulli.rst
213214
hypgeom.rst
214215
partitions.rst

‎examples/ode_basic.c‎

Lines changed: 69 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,69 @@
1+
#include <math.h>
2+
3+
#include "acb_types.h"
4+
#include "acb.h"
5+
#include "acb_mat.h"
6+
#include "acb_ode.h"
7+
#include "gr.h"
8+
#include "gr_ore_poly.h"
9+
10+
int main(void)
11+
{
12+
gr_ctx_t ZZ, Pol, Dop;
13+
gr_ptr dop;
14+
15+
int status = GR_SUCCESS;
16+
17+
slong prec = 50;
18+
19+
// Differential operator
20+
21+
// XXX fmpz_poly? fmpq_poly??
22+
// XXX std derivative
23+
gr_ctx_init_fmpz(ZZ);
24+
gr_ctx_init_gr_poly(Pol, ZZ);
25+
gr_ctx_init_gr_ore_poly(Dop, Pol, 0, ORE_ALGEBRA_EULER_DERIVATIVE);
26+
status |= gr_ctx_set_gen_name(Pol, "z");
27+
status |= gr_ctx_set_gen_name(Dop, "Tz");
28+
GR_TMP_INIT(dop, Dop);
29+
status |= gr_ore_poly_set_str(dop, "(z^2 + 1)*Tz^2 + (z^2 - 1)*Tz", Dop);
30+
// status |= gr_ore_poly_set_str(dop, "4*Tz^2 - 4*Tz - z^2 + 8*z - 11", Dop);
31+
// status |= gr_ore_poly_set_str(dop, "Tz^6 - 6*Tz^5 + 12*Tz^4 - 10*Tz^3 + 3*Tz^2 + z^2", Dop);
32+
// status |= gr_ore_poly_set_str(dop, "4*Tz^2 - 4*Tz - z^2 + 8*z - 11", Dop);
33+
GR_MUST_SUCCEED(status);
34+
35+
// flint_printf("dop = %{gr}\n", dop, Dop);
36+
37+
// Evaluation point
38+
39+
acb_t pt;
40+
acb_init(pt);
41+
// acb_set_d(pt, .5);
42+
acb_const_pi(pt, prec + 4);
43+
arb_one(acb_imagref(pt));
44+
acb_inv(pt, pt, prec + 4);
45+
46+
// Output matrix
47+
48+
acb_mat_t mat;
49+
slong dop_order = gr_ore_poly_length(dop, Dop) - 1;
50+
acb_mat_init(mat, dop_order, dop_order);
51+
52+
// Numerical solution
53+
54+
status |= acb_ode_fundamental_matrix(mat, dop, Dop, NULL, NULL, pt, 0, prec);
55+
GR_MUST_SUCCEED(status);
56+
57+
flint_printf("\n------\n");
58+
acb_mat_printd(mat, prec+4);
59+
flint_printf("\n%f\n", prec/log2(10.));
60+
61+
acb_clear(pt);
62+
acb_mat_clear(mat);
63+
GR_TMP_CLEAR(dop, Dop);
64+
65+
gr_ctx_clear(Dop);
66+
gr_ctx_clear(Pol);
67+
68+
flint_cleanup_master();
69+
}

‎examples/ode_fundamental_matrix_manual_exponents.c‎

Lines changed: 144 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,144 @@
1+
#include "acb_types.h"
2+
#include "acb.h"
3+
#include "acb_mat.h"
4+
#include "acb_ode.h"
5+
#include "gr.h"
6+
#include "gr_ore_poly.h"
7+
8+
// XXX also manual sing?
9+
10+
void
11+
fundamental_matrix(const char * dop_str,
12+
const acb_ode_exponents_struct * expos,
13+
double pt_d)
14+
{
15+
gr_ctx_t CC, Pol, Dop;
16+
gr_ptr dop;
17+
18+
slong prec = 30;
19+
20+
int status = GR_SUCCESS;
21+
22+
gr_ctx_init_complex_acb(CC, prec);
23+
gr_ctx_init_gr_poly(Pol, CC);
24+
gr_ctx_init_gr_ore_poly(Dop, Pol, 0, ORE_ALGEBRA_EULER_DERIVATIVE);
25+
26+
GR_TMP_INIT(dop, Dop);
27+
28+
status |= gr_ctx_set_gen_name(Pol, "z");
29+
status |= gr_ctx_set_gen_name(Dop, "Tz");
30+
status |= gr_ore_poly_set_str(dop, dop_str, Dop);
31+
32+
status |= gr_println(dop, Dop);
33+
flint_printf("\n");
34+
35+
GR_MUST_SUCCEED(status);
36+
37+
slong dop_order = gr_ore_poly_length(dop, Dop) - 1;
38+
39+
acb_mat_t mat;
40+
acb_mat_init(mat, dop_order, dop_order);
41+
42+
acb_t pt;
43+
acb_init(pt);
44+
acb_set_d(pt, pt_d);
45+
46+
GR_MUST_SUCCEED(acb_ode_fundamental_matrix(mat, dop, Dop, expos, NULL, pt, 0, prec));
47+
48+
flint_printf("%{acb_mat}\n\n", mat);
49+
flint_printf("--------\n\n");
50+
51+
acb_mat_clear(mat);
52+
GR_TMP_CLEAR(dop, Dop);
53+
gr_ctx_clear(Dop);
54+
gr_ctx_clear(Pol);
55+
gr_ctx_clear(CC);
56+
acb_clear(pt);
57+
}
58+
59+
60+
void
61+
apery(void)
62+
{
63+
acb_ode_shift_struct shift[1] = {{ .n = 0, .mult = 3 }};
64+
acb_ode_group_struct grp[1] = {{ .nshifts = 1, .shifts = shift }};
65+
acb_init(grp->leader);
66+
acb_zero(grp->leader);
67+
acb_ode_exponents_struct expos[1] = {{ .ngroups = 1, .groups = grp }};
68+
69+
fundamental_matrix(
70+
"(z^2 - 34*z + 1)*Tz^3 + (3*z^2 - 51*z)*Tz^2 + (3*z^2 - 27*z)*Tz + z^2 - 5*z",
71+
expos,
72+
0.015625);
73+
74+
acb_clear(grp->leader);
75+
}
76+
77+
78+
void
79+
multiple_shifts(void)
80+
{
81+
const char * dop = "Tz^6 - 6*Tz^5 + 12*Tz^4 - 10*Tz^3 + 3*Tz^2 + z^2";
82+
83+
acb_ode_shift_struct shift[3] = {
84+
{ .n = 0, .mult = 2 },
85+
{ .n = 1, .mult = 3 },
86+
{ .n = 3, .mult = 1 },
87+
};
88+
acb_ode_group_struct grp[1] = {{ .nshifts = 3, .shifts = shift }};
89+
acb_init(grp->leader);
90+
acb_zero(grp->leader);
91+
acb_ode_exponents_struct expos[1] = {{ .ngroups = 1, .groups = grp }};
92+
93+
fundamental_matrix(dop, expos, 2.);
94+
95+
acb_clear(grp->leader);
96+
}
97+
98+
99+
void
100+
whittaker(void)
101+
{
102+
slong prec = 30;
103+
104+
acb_t kappa, mu, half;
105+
acb_init(kappa);
106+
acb_init(mu);
107+
acb_init(half);
108+
109+
const char * dop = "4*Tz^2 - 4*Tz - z^2 + 8*z - 11";
110+
acb_set_si(kappa, 2);
111+
acb_set_si(mu, 3);
112+
acb_sqrt(mu, mu, prec);
113+
acb_set_d(half, .5);
114+
115+
acb_ode_shift_struct shift[1] = { { .n = 0, .mult = 1 } };
116+
acb_ode_group_struct grp[2] = {
117+
{ .nshifts = 1, .shifts = shift },
118+
{ .nshifts = 1, .shifts = shift },
119+
};
120+
acb_init(grp[0].leader);
121+
acb_sub(grp[0].leader, half, mu, prec);
122+
acb_init(grp[1].leader);
123+
acb_add(grp[1].leader, half, mu, prec);
124+
acb_ode_exponents_struct expos[1] = {{ .ngroups = 2, .groups = grp }};
125+
126+
fundamental_matrix(dop, expos, 2.);
127+
128+
acb_clear(grp[0].leader);
129+
acb_clear(grp[1].leader);
130+
acb_clear(kappa);
131+
acb_clear(mu);
132+
acb_clear(half);
133+
}
134+
135+
136+
int
137+
main(void)
138+
{
139+
apery();
140+
multiple_shifts();
141+
whittaker();
142+
143+
flint_cleanup_master();
144+
}

0 commit comments

Comments
 (0)