Re: Solving two-dimensional systems

[Bill Allombert <Bill.Allombert@math.u-bordeaux.fr>]

On Tue, May 26, 2015 at 08:24:48PM +0200, Bill Allombert wrote:
> On Tue, May 26, 2015 at 11:16:12AM -0400, Charles Greathouse wrote:
> > In computing https://oeis.org/A258042 I found myself solving this
> > two-dimensional system:
> > 
> > solve(a=1.7, 1.8, my(x=solve(y=1.8, 2,
> > y*(a+y)*log(a+y)-(a+y^2)*log(a+y^2))); log(a+x^2)/log(a+x)^2-1)
> > 
> > Is there any better way to do this than nested solve() calls?
> 
> First, you can simplify to
> 
> solve(a=1.7,1.8,my(x=solve(y=1.8,2,
>   y*(a+y)-(a+y^2)*log(a+y)));log(a+x^2)-log(a+x)^2)
> 
> Secondly, you could use a two dimensional Newton iteration:
> Set:
> F(a,x)=(x*(a+x)-(a+x^2)*log(a+x),log(a+x^2)-log(a+x)^2)
> 
> Compute the differential dF so that
> 
> F(a+h1,x+h2)=F(a,x)+dF(a,x)*[h1,h2]~+O(||(h1,h2)||^2)

Using diffop (I am lazy) I found:
F(x,a)=[x*(a+x)-(a+x^2)*log(a+x),log(a+x^2)-log(a+x)^2]~
dF(x,a)=my(L=log(x+a));[((-2*L+1)*x^2+(-2*L+3)*a*x+(a^2-a))/(x+a),((a-L)*x+(-L-1)*a)/(x+a);((-2*L+2)*x^2+2*a*x-2*L*a)/(x^3+a*x^2+a*x+a^2),(-2*L*x^2+x+(-2*L+1)*a)/(x^3+a*x^2+a*x+a^2)]
\p1000
V=[1.8,1.7]~;for(i=1,11,V-=matsolve(dF(V[1],V[2]),F(V[1],V[2])));V[2]
Cheers,
Bill.