wsolve.txt

这是Maple下使用的“吴方法”推理求解多项式方程组的源代码。

Leader := proc(dp, ord, _cls, _ord) 
         local lead, cls, order, L, l,  der, clsord, i, j, cls1, ord1,Diff_Var,Lvar; 
      global depend_relation; 
           options remember,system;     
           
#Step1 [Initialize] 

         Diff_Var:=depend_relation[2]; 
         Lvar:=nops(Diff_Var); 
         L:=indets(dp); l:=nops(L); 
         lead := 1; 
         cls := 0; 
         order := 0; 

#Step 2 [Algebraic mode] 

        if POLY_MODE='Algebraic' then 
         for i to nops(ord) do 
           if member(ord[i],L) then RETURN(ord[i]) fi 
         od; 
        fi; 

#Step2 [cls_ord_var case] 

         if RANKING_MODE = cls_ord_var then 
            for i from  1 to l do 
                der := L[i]; 

                cls1:=Class(der,ord); 
                ord1:=torder(der); 
                
                if cls1 > cls  or (cls1 = cls and ord1 > order) then 
                   lead := der; 
                   cls := cls1; 
                   order := ord1 
                else 
                   if cls1 = cls and ord1 = order and ord1 > 0 then 
                      for j from 1 to Lvar do 
                         if op(j, lead) < op(j, der) then 
                            lead := der; 
                            cls := cls1; 
                            order := ord1; 
                            break; 
                         else 
                            if op(j, lead) > op(j, der) then; 
                            break; 
                            fi; 
                         fi; 
                     od; 
                   fi; 
               fi; 
            od; 
         fi; 
          
#Step3 [ord_cls_var case] 
          
         if RANKING_MODE = ord_cls_var then 
            for i from  1 to l do 
                der := L[i]; 
                cls1:=Class(der,ord); 
                ord1:=torder(der); 
                if ord1 > order or (ord1 = order and cls1 > cls) then 
                   lead := der; 
                   cls := cls1; 
                   order := ord1 
                elif ord1 = order and cls1 = cls and ord1 > 0 then 
                     for j from 1 to Lvar do 
                         if op(j, lead) < op(j, der) then 
                            lead := der; 
                            cls := cls1; 
                            order := ord1; 
                            break; 
                         elif op(j, lead) > op(j, der) then; 
                            break; 
                         fi; 
                     od; 
               fi; 
            od; 
         fi; 
          

#STEP3 [Prepare for return] 

          
         if nargs > 2 then 
            _cls := cls; 
         fi; 
         if nargs > 3 then 
            _ord := order; 
         fi; 
          
         RETURN(lead); 
          
end: 



         


          
depend:=proc(l1,l2) 
global depend_relation; 
 depend_relation:=[l1,l2]; 
 1: 
 end: 
#============================================================ 
# 
#       {1,-1,0} <- dercmp(der1, der2) 
# 
#       input: der1, der2, two derivatives; 
#           
# 
#       output: 0 if der1 = der2 
#               1 if der1 < der2 
#              -1 if der1 > der2 
# 
#============================================================   
dercmp:=proc(der1, der2,ord) 
           local dp, lead; 
                options remember,system; 
            
           if der1=der2 then RETURN(0) fi; 
            
           dp := der1 + der2; 
           lead:=Leader(dp,ord); 
           if lead = der1 then 
              RETURN(-1); 
           else 
              RETURN(1); 
           fi; 
            
end: 

#if der1 is a derivative of der2 then return true else false
Is_derivative:=proc(der1, der2) 
              local dt1,dt2,i, l1,l2; 
                   options remember,system; 
              
              
               if type(der1,symbol) then dt1:=der1 else dt1:=op(0,der1) fi; 
               if type(der2,symbol) then dt2:=der2 else dt2:=op(0,der2) fi; 
            if dt1<>dt2 then RETURN(false) fi; 
            if der1=der2 then RETURN(true) fi; 
             
            if torder(der2)=0 then 
                 RETURN(true) 
            elif torder(der1)=0 then 
                   RETURN(false) ; 
           else 
                l1:=dOrder(der1,dt1); 
                l2:=dOrder(der2,dt1); 
                for i to nops(l1) do 
                     if op(i,l1) < op(i,l2) then RETURN(false) fi; 
                od; 
         fi; 
                           
           true; 
end: 
                                                  
#if der1 is a proper derivative of der2 then return true else false
Is_proper_derivative:=proc(der1, der2) 
              local dt1,dt2,i, l1,l2; 
                   options remember,system; 
              
              
               if type(der1,symbol) then dt1:=der1 else dt1:=op(0,der1) fi; 
               if type(der2,symbol) then dt2:=der2 else dt2:=op(0,der2) fi; 
            if dt1<>dt2 then RETURN(false) fi; 
            if der1=der2 then RETURN(false) fi; 
             
            if torder(der2)=0 then 
                 RETURN(true) 
            elif torder(der1)=0 then 
                   RETURN(false) ; 
           else 
                l1:=dOrder(der1,dt1); 
                l2:=dOrder(der2,dt1); 
                for i to nops(l1) do 
                     if op(i,l1) < op(i,l2) then RETURN(false) fi; 
                od; 
         fi; 
                           
           true; 
end: 

#return all the derivatives in dp of der
proper_derivatives:=proc(dp,der) 
    local i,idt,dets; 
         options remember,system; 
    idt:=indets(dp); 
    dets:={}; 
     
    for i in idt do 
      if Is_proper_derivative(i,der) then 
         dets:={i,op(dets)}; 
         
      fi; 
    od; 
    dets;      
end:      

       
     
     

#sometimes, Maple's GCD is NOT valid for parameters which has indexes such a[1],a[2]; 

NPremO := 

   proc(uu,vv,xx) 
   local r,v,dr,dv,lcv,rtv,g,lu,lv; 
    options remember,system; 
     if degree(vv,xx)=0 then RETURN(0) fi; 
     if type(vv/xx,integer) then coeff(uu,xx,0) 
       else 
           r := expand(uu); 
           dr := degree(r,xx); 
           v := expand(vv); 
           dv := degree(v,xx); 
            
             lcv:=lcoeff(v,xx); 
#            rtv:=sterm(v,xx); 
           rtv:=v-expand(lcoeff(v,xx)*xx^degree(v,xx)); 

           while dv <= dr and r <> 0 do 

#              g:=gcd(lcv,coeff(r,xx,dr)); 

             lu:=lcv; 

             lv:=coeff(r,xx,dr); 
                  
               r:=expand(lu*(r-expand(lcoeff(r,xx)*xx^degree(r,xx)))-lv*rtv*xx^(dr-dv)); 
               dr := degree(r,xx) 
           od; 
           r 
       fi 
   end: 
  
  
#######变量y有两种情形.y and y[1]所以要小心. 

  
ODPrem:= 

  proc(uu,vv,y) 
    local u,x,dv,du,d,t,var; 
    global depend_relation; 
    options remember,system;         
    
       t:=depend_relation[2][1]; 
       var:=depend_relation[1]; 
  
#        if dClass(uu,var)=0 then RETURN(uu) fi; 

          u:=uu; 
          if type(y,symbol) then x:=y else         x:=op(0,y)  fi;                 
          dv:=dOrder(vv,x); 
          du:=dOrder(u,x); 
          d:=du-dv; 
          if d<0 then RETURN(uu) fi; 
#the following program is for ordinary differential case                 
          while d>0 do 

            u:=NPremO(u,dDiff(vv,t,d), x[du]); 
            du:=dOrder(u,x); 
            d:=du-dv; 
          od; 
#         
          if dv =0 then NPremO(u,vv,x) else NPremO(u,vv,x[dv]) fi;                           

            
   end:     

    
NewPrem:= 

  proc(uu,vv,y) 
    local u,x,dv,du,d,t,var; 
    global depend_relation;
    options remember,system;      
    
    if POLY_MODE='Algebraic' then 
    
       NPrem(uu,vv,y); 
        
    elif POLY_MODE='Ordinary_Differential' then         
    
    ODPrem(uu,vv,y);                      
     
    elif POLY_MODE='Partial_Differential' then 
        
       PDPrem(uu,vv,y); 
        
    fi; 
            
   end:     



Headder:=proc(idts) 
local hder,i,j; 
     options remember,system; 
  hder:=op(1,idts);   
  
#   if type(hder,symbol) then var:=hder else var:=op(0,hder) fi; 
        
  for i in idts do 
  
    if torder(i) > torder(hder) then 
        hder:=i ;           
    elif torder(i)=torder(her) then 
       for j to nops(hder) do 
          if op(j,i) > op(j,hder) then hder:=i;break; 
          elif op(j,i) < op(j,hder) then break; 
          fi; 
       od; 
    fi; 
                   
   od; 
    
  hder: 
  
end: 


PDPrem:=proc(dp, dq, y) 
          local Diff_Indets,Diff_Var, var,dr, idt,l2,lder,da,l1,i,s; 
          options remember,system; 

#Step1 [Special case] 

          
    Diff_Indets:=depend_relation[1]; 
    Diff_Var:=depend_relation[2]; 
               
    if y=1 then RETURN(0) fi; 
     
    if type(y,symbol) then var:=y else var:=op(0,y) fi; 
     
    dr:=dp; 
     
       idt:=proper_derivatives(dr,y); 
       l2:=dOrder(y,var); 

#Step2 [Recursive base]         

       while  nops(idt)<>0 do 

          
          lder:=Headder(idt); 
          
          da := dq; 
                    

          
          l1:=dOrder(lder,var); 

          
          
          for i to nops(l1) do 
          
           if l2=0 then s:=op(i,l1) else s:=op(i,l1) -op(i,l2) fi; 
           
                  da := dDiff(da, Diff_Var[i],s); 
            
          od; 
          
          dr := NPremO(dr, da, lder); 
          
       idt:=proper_derivatives(dr,y);           
          
    od; 

         dr:=NPremO(dr,dq,y) ; 
          
#Step4 [Prepare for return] 

    dr;           
          
end:                 
          
  
    

# 
# set:  A set; 

cmb:=proc(set0) 
local p,q, ls,set1; 
     options remember,system; 
    ls:={}; 
    
    set1:=set0; 
    
 for p in set0 do 
    set1:=set1 minus {p}; 
    for q in set1  do 
        ls:={op(ls),[p,q]} 
    od; 
 od; 
 ls; 
end: 
  

######## modified Nov.6 ######## 

          
spolyset:=proc(as,ord) 
local res,bs,l,p,cls,i,poly; 
     options remember,system; 
    
   res:={}; 
   bs:={op(as)}; 
   l:=nops(bs); 
   if l <= 1 then RETURN(res) fi; 
    
   p:=op(1,bs); 
   cls:=Class(Leader(p,ord),ord); 
   bs:=bs minus {p}; 
        
   while nops(bs) <>0  do 
     for i in bs do 
        if Class(Leader(i,ord),ord) = cls then poly:=PD_spoly(i,p,ord); res:=res union {poly} fi; 
     od; 
      
      
   p:=op(1,bs); 
   cls:=Class(Leader(p,ord),ord); 
   bs:=bs minus {p};   
      
   od; 
  
   res minus {0}; 
end: 
    
          

      
            
    
    
PD_spoly:=proc(p,q,ord)     
local     leadp,leadq,var,l1,l2,ml,dp,dq,i,Diff_Var; 
global  depend_relation; 
options remember,system; 
    leadp:=Leader(p,ord); 
    leadq:=Leader(q,ord); 
     
    var:=op(0,leadp); 
    Diff_Var:=depend_relation[2]; 
     
    l1:=dOrder(leadp,var); 
    l2:=dOrder(leadq,var); 
     
    ml:=maxlist(l1,l2);