restart:
NOTE: please make sure that you are using a gfun version >= 3.66 to run these examples
gfun:-version();
# work around a 'feature' of gdev
`type/RootOf`:={'RootOf'(algebraic), 'RootOf'(algebraic, {range, constant, equation})};
interface(elisiondigitsthreshold=100,elisiondigitsbefore=10,elisiondigitsafter=10):
gfun
with(gfun);
deq[exp] := holexprtodiffeq(exp(z), y(z));
deq[Ai] := holexprtodiffeq(AiryAi(z), y(z));
`diffeq+diffeq`(deq[exp], deq[Ai], y(z));
rec := diffeqtorec(deq[Ai], y(z), u(n));
coeffproc := rectoproc(rec, u(n)):
coeffproc(50001);
cn := [1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796];
listtodiffeq(cn, y(z), ['ogf']);
listtoalgeq(cn, y(z), ['ogf']);
NumGfun
with(NumGfun);
Numerical Evaluation\302\240I:
A 'Less Common' Special Function
diffeq := collect({diff(y(z),z,z)
-(-2*z^5+4*z^3+z^4*a-2*z-a)/((z-1)*(z+1))^3*(diff(y(z),z))
-(-z^2*b+(-c-2*a)*z-d)/((z-1)^3*(z+1)^3)*y(z),
y(0)=1,D(y)(0)=0},diff,factor);
a, b, c, d := 1, 1/3, 1/2, 3;
evalf[51](HeunD(a, b, c, d, 1/3));
myHeunD := diffeqtoproc(diffeq, y(z)):
myHeunD(1/3, 50);
myHeunD(1/3, 400);
(More general code = less bugs!)
evalf(HeunD(a, b, c, d, -0.9));
myHeunD(-0.9, 9);
myHeunD(-0.9, 20);
evalf(HeunD(a, b, c, d, -0.99));
myHeunD(-0.99, 9);
myHeunD(-0.99, 100);
(No \342\200\234numerical unstability\342\200\235)
Numerical Evaluation\302\240II:
A Random Example
RandomTools:-SetState(state=42);
random_diffeq := proc(ord, d)
uses RandomTools; local myrat;
myrat := 'rational'('denominator'=60);
{ add(Generate('polynom'(myrat, z,
'degree'=d)) * diff(y(z), [z$k]),
k=0..ord),
seq((D@@k)(y)(0) = Generate(myrat),
k=0..ord-1) };
end proc:
rnd_diffeq := random_diffeq(4,3);
evaldiffeq(rnd_diffeq, y(z), 1/2, 30);
evaldiffeq(rnd_diffeq, y(z), (1+I)/3, 30);
evaldiffeq(rnd_diffeq, y(z), 1/3, 1000);
(Scales to high precisions!)
Analytic Continuation
diffeq := holexprtodiffeq(arctan(z), y(z));
diffeq := diffeqtohomdiffeq(diffeq, y(z));
evaldiffeq(diffeq, y(z), 1/2, 50);
barediffeq := op(1, diffeq);
evaldiffeq(barediffeq, y(z), (1+I)/2, 50);
transition_matrix(barediffeq, y(z), 1/2, 20);
evaldiffeq(diffeq, y(z), 5/4*(1+I), 50);
t := evaldiffeq(diffeq, y(z), [0,5/4*(1+I)], 50);
infolevel[gfun] := 2:
evaldiffeq(diffeq, y(z), [0,5/4*(1+I)], 50);
infolevel[gfun] := 0:
transition_matrix(diffeq, y(z), [0,I+1,2*I,I-1,0], 15);
Regular Singular Points
evaldiffeq(diffeq, y(z), [0,I], 20, ord=3);
Applications: Bessel functions, asymptotics, resummation, \342\200\234symbolic-numeric\342\200\235 differential algebra...
Regular Singular Points and Asymptotics
rec := { # Franel numbers, OEIS A000172
(n+2)^2*u(n+2)-(7*n^2+21*n+16)*u(n+1)-8*(n+1)^2*u(n),
u(1)=2, u(2)=10
};
rectoproc(rec, u(n), list)(10);
deq := collect(NumGfun:-utilities:-purerectodiffeq(rec, u(n), y(z)), diff, factor);
local_basis(deq, y(z), 0, 5);
local_basis(deq, y(z), 1/8, 5);
sol := analytic_continuation(deq, y(z), [0, 1/8], 10, ord=3);
sol := eval(sol, {_C[0]=0, _C[1]=1, ln=(x->ln(-x)+I*Pi)});
collect(simplify(equivalent(op([1,2], sol), z, n, 3)), n) assuming n::posint;
The 'Bit-Burst' Algorithm
diffeq := holexprtodiffeq(arctan(z), y(z));
my_e := convert(evalf[1000](exp(1)), rational, exact);
infolevel[gfun] := 2:
evaldiffeq(diffeq, y(z), [0,my_e], 1000);
infolevel[gfun] := 0:
TTdSMApJQVJUQUJMRV9TQVZFLzE4NDQ2ODg0MjE2MzIxNDA2MDA2WCwlKWFueXRoaW5nRzYiRiVbZ2whIiUhISEjJSIjIiMkIjYrKysrKysrKysrIiEjPyQiIiFGKiQiNUA7aCEzKzR3a2olRigkIjUrKysrKysrKyshKUYoRiU=TTdSMApJQVJUQUJMRV9TQVZFLzE4NDQ2ODg0MjE2Mzg3Njg3MjM4WCwlKWFueXRoaW5nRzYiRiVbZ2whIiUhISEjJSIjIiMkIjErKysrKysrNSEjOiQiIiFGKiQiMSR6KmVgRWZUSkYoRiZGJQ==