D\303\251mo NumGfun
S\303\251minaire Calin
15 mars 2011
restart: unassign('deq', 'rec'):
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);
NumGfun : \303\251valuation num\303\251rique
with(NumGfun);
Une fonction sp\303\251ciale pas si commune
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);
evalf(HeunD(a, b, c, d, -0.9));
myHeunD(-0.9, 9);
evalf(HeunD(a, b, c, d, -0.99));
myHeunD(-0.99, 400);
Un exemple \302\253\302\240au hasard\302\240\302\273
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);
Prolongement analytique
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), 1, 50);
t := evaldiffeq(diffeq, y(z), [0,1], 50);
infolevel[gfun] := 3:
evaldiffeq(barediffeq, y(z), [0,5/4*(1+I)], 20);
infolevel[gfun] := 0:
transition_matrix(diffeq, y(z), [0,I+1,2*I,I-1,0], 15);
Points singuliers r\303\251guliers
diffeq := collect(holexprtodiffeq(BesselI(sqrt(7),z), y(z)), diff, factor);
NumGfun:-transition_matrix(diffeq, y(z), [0,1/3], 20, 'forcepath');
ref := GAMMA(sqrt(7)+1)*exp(ln(2)*sqrt(7))*BesselI(sqrt(7), 1/3);
evalf[19](ref);
Asymptotique de suites r\303\251currentes
rec := {(n+2)^2*u(n+2)-(7*n^2+21*n+16)*u(n+1)-8*(n+1)^2*u(n), u(0)=2, u(1)=10};
deq := collect(NumGfun:-utilities:-purerectodiffeq(rec, u(n), y(z)), diff, factor);
sing := NumGfun:-utilities:-diffeq_infsing(deq, y(z), 'nonzero');
transition_matrix(deq, y(z), [0,2/5*1/8,3/5*1/8,1/8], 20, 'forcepath');
Compl\303\251ments
Bornes symboliques
R\303\251currences
rec := {(2*n+2)^2*u(n+1) = -(n+1)*u(n), u(0) = 1};
bound_rec(rec, u(n));
rec := {4*u(n+1)-13*u(n+2)+5*u(n+3), u(0) = 17/5, u(1) = 3, u(2) = 3/12, u(3) = 0, u(4) = 5};
bound_rec(rec, u(n));
\303\211quations diff\303\251rentielles
bound_diffeq({diff(y(z), z) = y(z), y(0) = 1}, y(z));
bound_diffeq(holexprtodiffeq(arctan(z), y(z)), y(z));
bound_diffeq(rnd_diffeq, y(z));
bound_diffeq(y(z) = (z^7+3*z^2+z+1)/((z-2)^3*(z^3-3)*(z-I)^2), y(z));
bound_diffeq_tail(holexprtodiffeq(arctan(z), y(z)), y(z), n);
diffeq := holexprtodiffeq(AiryAi(z), y(z));
dsolve(op(1,diffeq), y(z));
asympt(AiryBi(z), z,1);
bound_diffeq_tail(diffeq, y(z), 0);
Bornes num\303\251riques, contr\303\264le d'erreur
deq := holexprtodiffeq(arctan(z), y(z));
forget(NumGfun):
infolevel[gfun] := 3:
#infolevel[gfun] := 5:
evaldiffeq(deq, y(z), [0,11/10], 20);
infolevel[gfun] := 0:
R\303\251currences
Nombres de Catalan
rec := {C(n+1)=(4*n+2)/(n+2)*C(n),C(0)=1};
nth_term(rec, C(n), 50000);
Pi par la formule des Chudnovsky
A, B, C := 13591409, 545140134, 640320^3;
rec := {(A+B*n)*(3*n+3)*(3*n+2)*(3*n+1)*(n+1)^3*C * u(n+1) + (6*n+6)*(6*n+5)*(6*n+4)*(6*n+3)*(6*n+2)*(6*n+1)*(A+B*(n+1)) * u(n), u(0) = 12*A };
x := fnth_term(rec, u(n), 1000, 10001, 'series');
evalf[10001](sqrt(C)/x);
Points de grande taille
my_e := convert(evalf[1000](exp(1)), rational, exact);
infolevel[gfun] := 2:
evaldiffeq(diffeq, y(z), [0,my_e], 1000);
infolevel[gfun] := 0:
Approche des singularit\303\251s
infolevel[gfun] := 3:
evaldiffeq(barediffeq, y(z), 0.99*I, 20);
infolevel[gfun] := 0:
TTdSMApJNVJUQUJMRV9TQVZFLzc4NDU0OTEyWCwlKWFueXRoaW5nRzYiNiJbZ2whIiUhISEjJSIjIiMkIjYrKysrKysrKysrIiEjPyQiIiFGKyQiNUA7aCEzKzR3a2olRikkIjUrKysrKysrKyshKUYpRiY=TTdSMApJNVJUQUJMRV9TQVZFLzYzODkzNDAwWCwlKWFueXRoaW5nRzYiNiJbZ2whIiUhISEjJSIjIiMkIjErKysrKysrNSEjOiQiIiFGKyQiMSR6KmVgRWZUSkYpRidGJg==TTdSMApJNVJUQUJMRV9TQVZFLzg4NjIzNTA0WCwlKWFueXRoaW5nRzYiNiJbZ2whIiUhISEjJSIjIiMkIjdeSWUyek8mPl4kKnoiISM/JCE4JT4hUiopUWRoRiEqZVciRikkIjQnPmB0RiM+J2YyYkYpJCI1VC5kVCU0I2ZsJ1IlRilGJg==TTdSMApJNVJUQUJMRV9TQVZFLzU2NDAyNzIwWCwlKWFueXRoaW5nRzYiNiJbZ2whIiUhISEjJSIjIiMkIiIhRigkITZZbi9IN0I6KlI9QyEjPyQhNU1tOCd5JXBmX3ZPRiteJCQiNCE9SDE+amc5SFZGKyQiNi1IOkR6JFEwcWE2RitGJg==