Implicit: Taylor, Poiseux and Newton's diagram

Good day for all.

I wrote a few functions on the asymptotic expansion of implicit 
functions. I would be very interested to hear your opinion about my work:
1) how correctly my functions works;
2) whether they can be optimized and accelerated.

If these functions can be useful, then where and how I can put them for 
public use?

Evgeniy

/*==========================
X - list of variables
F - one expression (d=1) or a list of expressions of the variables X, d:length(F)
A - a fixed point, the value of X, F(A)=0
n - the order of the expansion

Let m:length(X). List X splits into two parts: X[1],...,X[m-d] - the independent variables,
X[m-d+1],...,X[m] - the dependent variables. If F at point A satisfies the implicit function theorem,
the imp_taylor(F,X,A,n) returns a single expansion or a list of expansions for X[m-d+1],...,X[m] in
powers of the differences (X[1]-A[1]), ..., (X[m-d]-A[m-d]) up to degree n inclusive.

Description of the method:
Let F: R^{m-d} * R^d -> R^d, F(A)=0. Let D is the Jacobian (matrix) of F at A. Then the only solution y(x) to the
equation F(x,y)=0 can be found by the method of successive approximations: y_{k+1}(x)=y_k(x)-D^{-1}F(x,y_k(x))

Examples:
imp_taylor(x-asin(y),[x,y],[%pi/6,1/2],5);
imp_taylor(x+s*sin(x+t)-x*exp(x-s),[t,s,x],[%pi,1,0],3);
imp_taylor([x[1]-log(y[1]),x[2]-exp(y[2])],[x[1],x[2],y[1],y[2]],[0,1,1,0],5);
===========================*/

imp_taylor(F,X,A,n):=block(
    [D,df,k],
	if listp(F) then block(
		[i,j,m:length(X),d:length(F)],
		if apply(ev,cons(F,map("=",X,A)))#makelist(0,d) then (
			print("Unable to compute taylor polynom: ",'at(F,X=A),"#",makelist(0,d)), 
			return(und)
		),
		genmatrix(lambda([i,j],diff(F[i],X[m-d+j])),d,d),
		D:apply(ev,cons(%%,map("=",X,A))),
		if is(determinant(D)=0) then (
			print("Unable to compute taylor polynom: determinant(D)=0"),
			return(und)
		),
		D:invert(D),
		F:taylor(apply(ev,cons(F,map("=",X,A+X))),X,0,n),
		df:makelist(0,d),
		for k:1 thru n do 
			df:df-transpose(D.map(expand,
taylor(apply(ev,cons(F,map("=",rest(X,m-d),df))),rest(X,-d),0,k) ))[1],
		rest(A,-d)+apply(ev,cons(df,map("=",X,X-A)))
	) else (
		if apply(ev,cons(F,map("=",X,A)))#0 then (
			print("Unable to compute taylor polynom: ",'at(F,X=A),"#0"), 
			return(und)
		),
		diff(F,last(X)),
		D:apply(ev,cons(%%,map("=",X,A))),
		if is(D=0) then (
			print("Unable to compute taylor polynom: ",'at('diff(F,last(X)),X=A),"=0"),
			return(und)
		),
		D:1/D,
		F:taylor(apply(ev,cons(F,map("=",X,A+X))),X,0,n),
		df:0,
		for k:1 thru n do 
			df:df-D*expand(taylor(ev(F,last(X):df),rest(X,-1),0,k)),
		last(A)+apply(ev,cons(df,map("=",X,X-A)))
	)
)$

/*==============================
X - list of 2 variables
F - an expression of the variables X
A - a fixed point, the value of X, F(A)=0
n - the order of the Poiseux expansion

Let F: R * R -> R satisfies the Taylor's theorem of n-th order. Let for simplicity X=[x,y]. Further suppose that:

at(F,[x=a,y=b])=0
at(diff(F,y,i),y=b)=0 for i=1,...,p-1
D: 1/p!* at(diff(F,y,p),[x=a,y=b])#0

Then the Poiseux's expansions of y(x) in a neighborhood of x=a can be found by the method of successive
approximations: y_{k+1}^p(x)=y_k^p(x)-D^{-1}F(x,y_k(x))

Examples:
imp_poiseux(x-y*exp(y),[x,y],[-%e^-1,-1],5);
imp_poiseux(y-2*sin(y/2)-sin(x)*(y^4+y^5)+cos(x)*(y^3-y^4),[x,y],[0,0],7);
===============================*/

imp_poiseux(F,X,A,n):=block(
    [df,F0,D,p:0,k,a,t],
    F0:ev(F,X[1]:A[1]),
    while (D:ev(F0,X[2]:A[2])/p!)=0 do (F0:diff(F0,X[2]), p:p+1),
    if is(p=0) then (
        print("Unable to compute: ",'at(F,[X[1]=A[1],X[2]=A[2]]),"#0"),
        return(und)
    ),
    F0:ratdisrep(taylor(X[2]^p-1/D*ev(F,X[1]:A[1]+X[1],X[2]:A[2]+X[2]),[X[1],0,n/p],[X[2],0,n])),
    if is(lopow(F0-ev(F0,X[2]:0),X[2])<p) then (
        print("Unable to compute: ",[X[1]=A[1],X[2]=A[2]]," is a singular point of ",F),
        return(und)
    ),
    df:0,
    for k:1 thru n do(
        ratdisrep(taylor(a*ev(F0,X[2]:df)^(1/p),X[1],0,k/p)),
        map(lambda([t],ratsubst(1,a^p,t)),%%),
        df:ratsimp(%%)
    ),
    print(a," in rootsof(",a^p-1,")"),
    A[2]+subst(X[1]=X[1]+A[1],expand(df))
)$

/*==============================
X - list of 2 variables
F - an irreducible polynomial of the variables X
A - a fixed point, the value of X, F(A)=0
n - the number of terms in the Poiseux expansion

Newton's diagram method for Poiseux's expansions of all branches of the algebraic function y(x)

Examples:
imp_newton(x^2+y^2-1,[x,y],[1/2,sqrt(3)/2],5)$
imp_newton(2*x^7-x^8-0*x^3*y+(4*x^2+x^3)*y^2+(x^3-x^4)*y^3-4*x*y^4+7*x^5*y^5+(1-x^2)*y^6+5*x^6*y^7+x^3*y^8,[x,y],[0,0],3)$
imp_newton(ratdisrep(taylor(x-asin(y),y,0,9)),[x,y],[0,0],5)$
===============================*/

imp_newton(F,X,A,n):=block(
    [F0,L,P,h,i],
    local(S),

    S(L,P):=block(
        [i,j,a,b,q,k,m,t,p:[],c:[],Q:[],D:[],d,P0,l],
        for i:1 thru h+1 do 
            if P[i]#0 then (a:lopow(P[i],X[1]), p:endcons([a,i-1],p), c:endcons(coeff(P[i],X[1],a),c)),
        i:1,
        while i<length(p) do (
            q:[], a:minf, 
            for j:length(p) step -1 thru i+1 do (
                b:(p[i][1]-p[j][1])/(p[j][2]-p[i][2]),
                if b=a then q:cons(j,q) elseif b>a then (a:b, q:[j])
            ),
            q:cons(i,q),
            Q:endcons(q,Q),
            D:endcons(a,D),
            i:last(q)
        ),
        if length(L)=0 then d:minf else d:last(L)[1],
        for i:1 thru length(Q) while (D[i]>d and D[i]#0) do
            for b in map(rhs,solve( sum(c[Q[i][j]]*m^(p[Q[i][j]][2]-p[Q[i][1]][2]),j,1,length(Q[i])) ,m)) do (
                P0:makelist(expand(sum(k!/(j!*(k-j)!)*b^(k-j)*X[1]^(D[i]*(k-j))*P[k+1],k,j,h)),j,0,h),
                l:endcons([D[i],b],L),
                if length(l)<n then S(l,P0) else print(A[2]+lsum(t[2]*(x-A[1])^t[1],t,l))
            )
    ),

    if is(ev(F,X[1]:A[1],X[2]:A[2])#0) then (
        print("Unable to compute: ",'at(F,[X[1]=A[1],X[2]=A[2]]),"#0"),
        return(und)
    ),

    F0:ev(F,X[1]:A[1]+X[1],X[2]:A[2]+X[2],expand),
    h:hipow(F0,X[2]),
    P:makelist(coeff(F0,X[2],i),i,0,h),

    S([],P)
)$

_______________________________________________
Maxima mailing list
Maxima <at> math.utexas.edu
http://www.math.utexas.edu/mailman/listinfo/maxima