# This is a program for real root isolation of multivariate real nonlinear systems # Before use our program, one need to load the Maple package intpakX for interval computation # We provide two commands: # 1: [R,SR]=rootfinding_cotrans(sys,var,B0,err,k) # 2: [R,SR]=rootfindingbox(sys,var,B0,err,k) # The first one is based on coordinate transformation and then we need only to find real roots in subboxes of [-1,1]*[-1,1] # The second one is that finding real roots in the given box B0 ####################### # Parameters: # Input: # sys: A multivariate real nonlinear systems (f1,f2,...,fn) # var: list of variables # B0: An initial box # err: A termination precision and err>0 # k: (optional, default=0) An integer for spliting the initial box into 2^k parts before computation # Output: # R: The isolating box set of the given systen # SR: The suspected root box set of the given systen ###################### # Examples: # libname:="D:\\intpakX", libname; # read("D:\\rootfinding_analyticfun.txt"); #f1:=exp(x+y)+x^2-3;;f2:=cos(y)-x^2; #A:=rootfinding_cotrans([f1,f2],[x,y],[[-10,10],[-10,10]],10^(-4)); # [[[[0.750000000000000, 1.00000000000000], # [-0.500000000000000, -0.250000000000000]], [ # [-0.500000000000000, -0.375000000000000], # [1.39130434782609, 1.50943396226415]]], []] # R:=A[1];SR:=A[2]; with(RootFinding): with(ListTools): with(VectorCalculus): with(LinearAlgebra): with(intpakX): with(combinat): Digits:=15: rootfinding_cotrans:=proc(sys,var,B0,err) local uB,n,Q,T,i,j,B,R,SR,t,L,temp,F,k; if nargs<5 then k:=0; else k:=args[5]; fi; if nops(sys)<>2 or nops(var)<>2 then ERROR("The input system must be a one with two functions and two variables!"); fi; R:=[];SR:=[]; L:=allbinary(n); uB:=unitbox(B0); for i from 1 to 4 do temp:=transformsys(sys,var,uB,L[i]); F:=temp[1];B:=temp[2]; for j from 1 to nops(B) do if k=0 then T:=rootfindingbox(F,var,B[j],err); else T:=rootfindingbox(F,var,B[j],err,k); fi; Q:=transformbox(T[1],L[i]); R:=[op(R),op(Q)]; Q:=transformbox(T[2],L[i]); SR:=[op(SR),op(Q)]; od; od; T:=SR;SR:=[]; for i from 1 to nops(T) do Q:=rootfindingbox(sys,var,T[i],err); R:=[op(R),op(Q[1])]; SR:=[op(SR),op(Q[2])]; od; return [R,SR]; end: transformbox:=proc(B,L) local n,Q,T,i,j,t; n:=nops(L);Q:=[]; for i from 1 to nops(B) do; T:=B[i]; for j from 1 to nops(L) do if L[j]=1 then t:=[1/T[j][2],1/T[j][1]]; T[j]:=t; fi; od; Q:=[op(Q),T]; od; return Q; end: unitbox:=proc(B) local n,i,V,Q,T; n:=nops(B);V:=[]; for i from 1 to n do T:=B[i]; if T[1]<-1 and T[2]>1 then Q:=[[-1,1/T[1]],[-1,1],[1/T[2],1]]; elif T[1]<-1 then if T[2]<=-1 then Q:=[[1/T[2],1/T[1]],[],[]]; else Q:=[[-1,1/T[1]],[-1,T[2]],[]]; fi; elif T[2]>1 then if T[1]>=1 then Q:=[[],[],[1/T[2],1/T[1]]]; else Q:=[[],[T[1],1],[1/T[2],1]]; fi; else Q:=[[],T,[]]; fi; V:=[op(V),Q]; od; return V; end: transformsys:=proc(sys,var,B,L) local n,F,T,Q,t,i,j,k,a; n:=nops(sys);F:=sys;Q:=[];k:=1; if L[1]=0 then if B[1][2]=[] then k:=0; else Q:=[[B[1][2]]]; fi; fi; if L[1]=1 then F:=subs(var[1]=1/var[1],F); if B[1][1]=[] and B[1][3]=[] then k:=0 else if B[1][1]<>[] then Q:=[[B[1][1]]]; fi; if B[1][3]<>[] then Q:=[op(Q),[B[1][3]]]; fi; fi; fi; i:=2; while i<=n and k=1 do if L[i]=0 then if B[i][2]=[] then k:=0; else T:=[]; for j from 1 to nops(Q) do a:=[op(Q[j]),B[i][2]]; T:=[op(T),a]; od; Q:=T; fi; fi; if L[i]=1 then F:=subs(var[i]=1/var[i],F); if B[i][1]=[] and B[i][3]=[] then k:=0 else T:=[]; if B[i][1]<>[] then for j from 1 to nops(Q) do a:=[op(Q[j]),B[i][1]]; T:=[op(T),a]; od; fi; if B[i][3]<>[] then for j from 1 to nops(Q) do a:=[op(Q[j]),B[i][3]]; T:=[op(T),a]; od; fi; Q:=T; fi; fi; i:=i+1 od; F:=simplify(F); F:=numer(F); if k=0 then Q:=[]; fi; Q:=evalf(Q); return [F,Q]; end: transformsys_old:=proc(sys,var,N,L) local n,F,B,T,Q,t,i,j; n:=nops(sys);F:=sys;B:=[[seq([-1.,1.],i=1..n)]]; for i from 1 to n do if L[i]=1 then F:=subs(var[i]=1/var[i],F); Q:=[]; for j from 1 to nops(B) do T:=B[j]; t:=evalf([-1,-1/N]);T[i]:=t; Q:=[op(Q),T]; t:=evalf([1/N,1]);T[i]:=t; Q:=[op(Q),T]; od; B:=Q; fi; od; F:=simplify(F); return [F,B]; end: allbinary:=proc(n) local i,j,T,Q; T:=[[0],[1]]; for i from 2 to n do Q:=[]; for j from 1 to nops(T) do Q:=[op(Q),[op(T[j]),0]]; Q:=[op(Q),[op(T[j]),1]]; od; T:=Q; od; return T; end: rootfindingbox:=proc(sys,var,B0,err) local BS,T,R,SR,n,a1,a2,a3,w,B,F,m,sys_horner,F_horner; if nargs<5 then BS:=evalf([B0]); else BS:=splitB(evalf(B0),args[5]); fi; R:=[];SR:=[]; n:=nops(BS);m:=nops(BS[1]); sys_horner:=convert(sys,horner,var); while n>0 do B:=BS[1]; BS:=subsop(1=NULL,BS); a1:=Exclusion(sys_horner,var,B); if a1=0 then F:=preconditioner(sys,var,B); #F_horner:=convert(F,horner,var); #a2:=IsOMsys(F_horner,var,B); a2:=IsOMsys(F,var,B); if a2=1 then a3:=Existence(F,var,B); if a3=1 then R:=[op(R),B]; elif a3=1/2 then SR:=[op(SR),B]; fi; else w:=widthB(B,m); if w>err then T:=splitB(B,2); BS:=[op(BS),op(T)]; else SR:=[op(SR),B]; fi; fi; fi; n:=nops(BS); od; return [R,SR]; end: middleB:=proc(B,n) local p,i; p:=[]; for i from 1 to n do p:=[op(p),evalf((B[i][1]+B[i][2])/2)]; od; return p; end: widthB:=proc(B,n) local p,i; p:=[]; for i from 1 to n do p:=[op(p),B[i][2]-B[i][1]]; od; p:=sort(p); return p[n]; end: mysubs:=proc(f,var,p) local A,a,i; A:=[]; for i from 1 to nops(p) do A:=[op(A),var[i]=p[i]]; od; a:=subs(op(A),f); return evalf(a); end: intervalfun:=proc(f,var,B) local F,A; if type(f,numeric) then A:=evalf([f,f]); else F:=inapply(f,op(var)); A:=F(op(B)); fi; return A; end: intervalfun_mean:=proc(f,var,B) local n,p,fp,A,a,i; n:=nops(B);p:=middleB(B,n); fp:=mysubs(f,var,p); a:=intervalfun(diff(f,var[1]),var,B) &* (B[1] &- p[1]); A:=fp &+ a; for i from 2 to n do a:=intervalfun(diff(f,var[i]),var,B) &* (B[i] &- p[i]); A:=A &+ a; od; return A; end: Exclusion:=proc(F,var,B) local ind,i,A; ind:=0;i:=1; while i<=nops(F) and ind=0 do A:=intervalfun(F[i],var,B); if A[1]>0 or A[2]<0 then ind:=1; fi; i:=i+1; od; return ind; end: generateOMmatrix:=proc(n) local i,U,k; k:=0; while k=0 do U:=RandomMatrix(n,generator=-1. .. 1. ); for i from 1 to n do U[[i],[i]]:=sign(U[i][i])*n; od; k:=IsOMmatrix(U,n); od; return U; end: Isnumericmatrix:=proc(M,n) local i,j,t; t:=1;i:=1;j:=1; while i<=n and t=1 do while j<=n and t=1 do if type(M[i][j],numeric) then else t:=0; fi; j:=j+1; od; i:=i+1; od; return t; end: preconditioner:=proc(sys0,var,B) local p,Jp,A,sys,U,n,a; n:=nops(B); p:=middleB(B,n); Jp:=Jacobian(sys0 ,var = p); if Rank(Jp)=n then if n=2 then U:=Matrix(2,2,[[2,1],[1,-2]]); elif n=3 then U:=Matrix(3,3,[[3,1,1],[1,-3,1],[1,1,3]]); elif n=4 then U:=Matrix(4,4,[[4,1,1,-1],[1,-4,1,1],[1,1,4,1],[1,1,1,4]]); elif n=5 then U:=Matrix(5,5,[[5,1,1,-1/2,1],[1,-5,1,-1,-1],[-1,1,5,1,-1],[1,1,1,5,-1],[1,1,1,1,5]]); else U:=generateOMmatrix(n); fi; A:=U.MatrixInverse(Jp); a:=Isnumericmatrix(A,n); if a=0 then A:=Matrix(n, shape = identity); fi; else A:=Matrix(n, shape = identity); fi; sys:=A.; sys:=convert(sys,list); return sys; end: IsOMmatrix:=proc(M,n) local ind,i,j,W,T; ind:=1;i:=1; while i<=n and ind=1 do j:=1; while j<=n and ind=1 do if M[i][j]=0 then ind:=0; fi; j:=j+1; od; i:=i+1; od; i:=n; while i>=2 and ind=1 do T:=choose(n,i); j:=1; while j<=nops(T) and ind=1 do W:=Determinant(M[[1..i],T[j]]); if W=0 then ind:=0; fi; j:=j+1; od; i:=i-1; od; return ind; end: IsOMsys:=proc(F,var,B) local ind,i,j,S,T,B0,f,svar,W,M,J,n; ind:=1;i:=1;J:=Jacobian(F,var);n:=nops(F); M:=Matrix(n,n); while i<=n and ind=1 do j:=1; while j<=n and ind=1 do M[[i],[j]]:=intervalfun(J[i][j],var,B); if M[i][j][1]<=0 and M[i][j][2]>=0 then ind:=0; fi; j:=j+1; od; i:=i+1; od; i:=n; while i>=2 and ind=1 do S:=Matrix(i,i,symbol=s); f:=Determinant(S); svar:=convert(S,list); T:=choose(n,i); j:=1; while j<=nops(T) and ind=1 do B0:=convert(M[[1..i],T[j]],list); W:=intervalfun(f,svar,B0); if W[1]<=0 and W[2]>=0 then ind:=0; fi; j:=j+1; od; i:=i-1; od; return ind; end: splitB:=proc(B,z) local n,T,m,Q,temp,i,j,a,k,w; n:=nops(B);T:=[B]; for i from 1 to n do m:=nops(T);w:=(B[i][2]-B[i][1])/z;Q:=[]; for j from 1 to m do for k from 0 to z-1 do; a:=[B[i][1]+k*w,B[i][1]+(k+1)*w]; temp:=subsop(i=a,T[j]); Q:=[op(Q),temp]; od; od; T:=Q; od; return T; end: vertexB:=proc(B) local p,n,P,i,j; n:=nops(B);p:=[[B[1][1]],[B[1][2]]]; for i from 2 to n do P:=[]; for j from 1 to nops(p) do P:=[op(P),[op(p[j]),B[i][1]]]; P:=[op(P),[op(p[j]),B[i][2]]]; od; p:=P; od; return p; end: mysubs2:=proc(f,var,p) local n,q,ind,i; n:=nops(p); q:=mysubs(f,var,p[1]); if q=0 then ind:=0; else ind:=sign(q); fi; i:=2; while i<=n and ind<>0 do q:=mysubs(f,var,p[i]); if q=0 then ind:=0; elif sign(q)<>ind then ind:=0; fi; i:=i+1; od; return ind; end: enlargeB:=proc(B0,B) local n,w1,w2,A,i; n:=nops(B);A:=[]; for i from 1 to n do w1:=(B0[i][1]-B[i][1])/2; w2:=(B[i][2]-B0[i][2])/2; A:=[op(A),[B0[i][1]-w1,B0[i][2]+w2]]; od; return A; end: Existence:=proc(F,var,B) local n,i,T,G,B0,var0,vB,s1,s2,ind,A,w,setB,bigB,k; n:=nops(B); if n=2 then A:=Existence_2d(F,var,B); ind:=A[1]; else i:=n;T:=[];k:=0; #### find two systems in faces #### while i>=1 and nops(T)<2 do G:=subs(var[i]=B[i][1],F); B0:=subsop(i=NULL,B); var0:=subsop(i=NULL,var); A:=Existence(G[1..n-1],var0,B0); if A=1 then T:=[op(T),[G,var0,B0]]; elif A=1/2 then k:=1; fi; G:=subs(var[i]=B[i][2],F); A:=Existence(G[1..n-1],var0,B0); if A=1 then T:=[op(T),[G,var0,B0]]; elif A=1/2 then k:=1; fi; i:=i-1; od; ##### compute teh evaluation of fn in two faces ###### if nops(T)=2 then G:=T[1][1];var0:=T[1][2];B0:=T[1][3]; w:=widthB(B0,n-1);setB:=[B0]; vB:=vertexB(B0); s1:=mysubs2(G[n],var0,vB); while s1=0 and w>10^(-6) do i:=1;setB:=splitB(B0,2);A:=0; while i<=nops(setB) do A:=Existence(G[1..n-1],var0,setB[i]); if A=1 then B0:=setB[i]; i:=nops(setB)+1; elif A=1/2 then bigB:=enlargeB(setB[i],B0); A:=Existence(G[1..n-1],var0,bigB); if A=1 then B0:=bigB; i:=nops(setB)+1; fi; fi; i:=i+1; od; if A=1 then vB:=vertexB(B0); s1:=mysubs2(G[n],var0,vB); if s1=0 then w:=w/2; fi; else w:=10^(-7); fi; od; G:=T[2][1];var0:=T[2][2];B0:=T[2][3]; w:=widthB(B0,n-1);setB:=[B0]; vB:=vertexB(B0); s2:=mysubs2(G[n],var0,vB); while s2=0 and w>10^(-6) do i:=1;setB:=splitB(B0,2);A:=0; while i<=nops(setB) do A:=Existence(G[1..n-1],var0,setB[i]); if A=1 then B0:=setB[i]; i:=nops(setB)+1; elif A=1/2 then bigB:=enlargeB(setB[i],B0); A:=Existence(G[1..n-1],var0,bigB); if A=1 then B0:=bigB; i:=nops(setB)+1; fi; fi; i:=i+1; od; if A=1 then vB:=vertexB(B0); s2:=mysubs2(G[n],var0,vB); if s2=0 then w:=w/2; fi; else w:=10^(-7); fi; od; if s1*s2<0 then ind:=1; elif s1*s2>0 then ind:=0; else ind:=1/2; fi; elif k=0 then ind:=0; else ind:=1/2; fi; fi; return ind; end: Findindex:=proc(A,a) local n,i,t; n:=nops(A);i:=1;t:=0; while i<=n and t=0 do if A[i]=a then t:=i; fi; i:=i+1; od; return t; end: whichcase:=proc(f,g) local sf,sg,t,ind; sf:=sort(f); sg:=sort(g); if sf[1]>0 or sf[4]<0 or sg[1]>0 or sg[4]<0 then ind:=0; elif (sf[1]=0 and sf[2]>0) or (sf[4]=0 and sf[3]<0) then t:=Findindex(f,0); if g[t]=0 then ind:=1; else ind:=0; fi; elif (sg[1]=0 and sg[2]>0) or (sg[4]=0 and sg[3]<0) then t:=Findindex(g,0); if f[t]=0 then ind:=1; else ind:=0; fi; else ind:=2; fi; return [ind,t]; end: bisection:=proc(f,g,var,L0,F0,G0) local F,G,s,L,w,fp,gp,p; F:=F0;G:=G0;s:=0;L:=L0; if G[1]*G[2]>0 then s:=G[1]; else w:=abs(L[2]-L[1]); while w>10^(-10) and s=0 do p:=(L[1]+L[2])/2; fp:=evalf(subs(var=p,f)); gp:=evalf(subs(var=p,g)); w:=abs(L[2]-L[1]); if F[1]*fp<=0 and G[1]*gp<=0 then L[2]:=p; elif F[2]*fp<=0 and G[2]*gp<=0 then L[1]:=p; else s:=gp; fi; od; fi; return s; end: Existence_2d:=proc(sys,var,B) local f,g,uf,ug,p1,p2,p3,p4,f1,f2,f3,f4,g1,g2,g3,g4,s,c,t,ind,P,G,F,L,S; f:=sys[1];g:=sys[2]; p1:=[B[1][1],B[2][2]]; p2:=[B[1][2],B[2][2]]; p3:=[B[1][2],B[2][1]]; p4:=[B[1][1],B[2][1]]; P:=[p1,p2,p3,p4]; f1:=evalf(subs(var[1]=p1[1],var[2]=p1[2],f)); f2:=evalf(subs(var[1]=p2[1],var[2]=p2[2],f)); f3:=evalf(subs(var[1]=p3[1],var[2]=p3[2],f)); f4:=evalf(subs(var[1]=p4[1],var[2]=p4[2],f)); g1:=evalf(subs(var[1]=p1[1],var[2]=p1[2],g)); g2:=evalf(subs(var[1]=p2[1],var[2]=p2[2],g)); g3:=evalf(subs(var[1]=p3[1],var[2]=p3[2],g)); g4:=evalf(subs(var[1]=p4[1],var[2]=p4[2],g)); F:=[f1,f2,f3,f4];G:=[g1,g2,g3,g4]; c:=whichcase([f1,f2,f3,f4],[g1,g2,g3,g4]); ################################################ if c[1]=0 then ind:=0; elif c[1]=1 then ind:=1/2; t:=c[2]; else S:=[]; if f1=0 then f1:=1; fi; if f2=0 then f2:=1; fi; if f3=0 then f3:=1; fi; if f4=0 then f4:=1; fi; if f1*f2<0 then L:=B[1]; uf:=subs(var[2]=B[2][2],f); ug:=subs(var[2]=B[2][2],g); s:=bisection(uf,ug,var[1],L,F[1..2],G[1..2]); S:=[op(S),s]; if s=0 then t:=12; fi; fi; if f2*f3<0 then L:=[B[2][2],B[2][1]]; uf:=subs(var[1]=B[1][2],f); ug:=subs(var[1]=B[1][2],g); s:=bisection(uf,ug,var[2],L,F[2..3],G[2..3]); S:=[op(S),s]; if s=0 then t:=23; fi; fi; if f3*f4<0 then L:=[B[1][2],B[1][1]]; uf:=subs(var[2]=B[2][1],f); ug:=subs(var[2]=B[2][1],g); s:=bisection(uf,ug,var[1],L,F[3..4],G[3..4]); S:=[op(S),s]; if s=0 then t:=34; fi; fi; if f4*f1<0 then L:=B[2]; uf:=subs(var[1]=B[1][1],f); ug:=subs(var[1]=B[1][1],g); s:=bisection(uf,ug,var[2],L,[F[4],F[1]],[G[4],G[1]]); S:=[op(S),s]; if s=0 then t:=41; fi; fi; if nops(S)>1 then if S[1]*S[2]<0 then ind:=1; elif S[1]*S[2]>0 then ind:=0; else ind:=1/2; fi; else ind:=1/2; t:=0; fi; fi; return [ind,t]; end: