wsolve.txt

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

#      	  print(nn);
       if Class(i[nops(i)],ord)=1 and POLY_MODE='Algebraic' then Is:=Inits1(i,ord) else Is := Inits(i,ord) fi; 
       if POLY_MODE='Partial_Differential' or POLY_MODE='Ordinary_Differential' then Is := Is union Septs(i,ord) fi; 
       Ninits := Is minus {op(nzero)}; 
#        print("checknum=",checknum, "Ninits=",Ninits); 
       for j to nops(Ninits) do 
           ps1 := ({op(ps)} union {op(i)}) union {Ninits[j]}; 
           if j = 1 then nzero1 := {op(nzero)} 
           else nzero1 := nzero1 union {Ninits[j-1]} 
           fi; 
           cset := cset union Css(ps1,ord,nzero1,T) 
       od 
   od; 
   cset 
end: 

realfac:=proc(ps,ord) 
local  s,res,i; 
     options remember,system; 
      s:=Produce(ps,ord,{}); 
      res:={}; 
      for i in s do 
         res:=res union i; 
      od; 
      res; 
end:               
  
Degenerate:=proc(ds,as,ord) 
local i,r; 
     options remember,system; 
       for i in ds do 
         r:=Premas(i,as,ord,'std_asc'); 
         if r =0 then return 1 fi; 
       od; 
         0; 
end: 
          
###  POLY_MODE:= one of ['Algebraic','Ordinary_Differential','Partial_Differential']
### depend_relation  is like : [[x,y],[t]];          
#### T:= one of ['std_asc','reg_asc','wu_asc','gao_asc','tri_asc']          
CSS:=proc(ps,ord,nzero,T) 
global factor_time,prem_time, t_time,time_bs,time_rs,time_f; 
local asc,i,pset,cset,nonzero,order,result,wm; 
options remember,system; 
   if nargs < 2 then order:=[op(indets(ps))] else order:=[op(ord)] fi; 
   if nargs < 3 then nonzero:={} else nonzero:=realfac({op(nzero)},ord) fi; 
   if nargs < 4 then asc:='std_asc' else asc:=T fi; 
   if asc='reg_asc' then return Dec({op(ps)},order,nonzero,'std_asc');fi; 
   cset:={}; 
   factor_time:=0; 
   prem_time:=0; 
   time_bs:=0; 
   time_rs:=0; 
   time_f:=0; 
    
   t_time:=time(); 
   pset:=Nrs(ps,order,nonzero); 

   
   if asc='reg_asc' then 
        for i to nops(pset) do cset:=cset union Dec(pset[i],order,nonzero,'std_asc') od;
   else         
   	for i to nops(pset) do cset:=cset union Css(pset[i],order,nonzero,asc) od; 
   fi;
   result:=[]; 
   for i in cset do if Degenerate(nonzero,i,order)<>1 then result:=[op(result),i]; fi od; 
       lprint("Timing","Total", time()-t_time, "Factoring", factor_time,"Prem",prem_time); 
       lprint("BasicSet",time_bs,"RS",time_rs,"factor",time_f); 
   result; 

end: 

charset:=proc(ps,ord,nzero,T) 
global factor_time,prem_time, t_time,time_bs,time_rs,time_f; 
local asc,nonzero,order,result,cm; 
options remember,system; 
   if nargs < 2 then order:=[op(indets(ps))] else order:=[op(ord)] fi; 
   if nargs < 3 then nonzero:={} else nonzero:=realfac({op(nzero)},ord) fi; 
   if nargs < 4 then asc:='std_asc' else asc:=T fi; 

   prem_time:=0; 
   time_bs:=0; 
    
  
        
   t_time:=time(); 

   result:=charset_b(map(primpart,{op(ps)},order),order,nonzero,asc); 

       lprint("Timing","Total", time()-t_time, "BasicSet",time_bs,"Prem",prem_time); 

   result; 

end: 

Css1:=proc(ps,ord,test) 
local  cset,i,j,septs,cset1,vv; 
     options remember,system; 
       cset:=CSS(ps,ord,test,'std_asc'); 
#       septs:={}; 
       cset1:=cset; 
       for i in cset1 do 
               septs:=Septs(i,ord) minus Inits(i,ord); 
               if septs<>{} then

                       for j in septs do 
                       vv:=  {op(Css1(ps union {op(i)} union {j},ord,test))};                      
                       cset:={op(cset)} union vv ;
                       od 
               fi 
       od; 
   [op(cset)] 
end: 

Gsolve:=proc(ps,ord,test) 
    gsolve(ps,ord,test) 
end: 

e_val:=proc(ps,ord,test) 
    Css1({op(ps)},ord,test) 
end: 

 

wsolve:=proc() 
local rt; 
global POLY_MODE; 
POLY_MODE:='Algebraic'; 
 rt:=CSS(args); 
 rt; 
end: 

od_wsolve:=proc() 
local rt; 
global POLY_MODE,depend_relation; 
POLY_MODE:='Ordinary_Differential'; 
 if type(depend_relation,'list') and nops(depend_relation) >1 and 
    type(depend_relation[2],'list') and nops(depend_relation[2])=1 then 
     rt:=CSS(args); 
  else 
     lprint("Please set depend_relation:=[[x,y,...],[t]] first, if x,y,... depend on t"); 
     rt:=0; 
  fi; 
 rt; 
end: 

pd_wsolve:=proc() 
local rt; 
global POLY_MODE,depend_relation; 
POLY_MODE:='Partial_Differential'; 
 if type(depend_relation,'list') and nops(depend_relation) >1 and 
    type(depend_relation[2],'list') and nops(depend_relation[2])>1 then 
     rt:=CSS(args); 
  else 
     lprint("Please set depend_relation:=[[x1,x2,...],[t1,t2,...]] first, if x1,x2,... depend on t1,t2,..."); 
     rt:=0; 
  fi; 
 rt; 
end: 

INDEX:=0: 
__testp__:=0:
remfac:={}: 
prem_time:=0: 
factor_time:=0: 
time_bs:=0: 
time_rs:=0: 
time_f:=0: 
POLY_MODE:=Algebraic: 
RANKING_MODE:='cls_ord_var': 

# Test Program for Differentiation 

#---------------------------------------- 
# Global dependence 
#---------------------------------------- 

######################################################### 
#11111111111111111111111111111111111111111111111111111111 
# 
#Part IV: The basic functions for differential poly 
#         
######################################################### 

#Diff_Var:=[u,v,w]; # Declare variables; 

#Diff_Indets:=[X1,X2,X3]; # Declare Diff ndeterminates; 

#Lvar:=nops(Diff_Var); 

#---------------------------------------------------------- 
# main function 
#---------------------------------------------------------- 

dDiff:=proc(dp,var,n) 
   local df,i; 
     options remember,system;    
   df:=dp; 
   if nargs=2 then RETURN (DF(dp,var)) fi; 
   for i to n do 
      df:=DF(df,var); 
   od; 
   df; 
end: 
          
#========================================================== 
#         
#          dq  <- DF(dp, var) 
# 
#          (differetiating a diff polynomial) 
# 
# 
#          Input: dp, a diff polynomial; 
#                 var,  the variable w.r.t which to be differentiate   
# 
#          Output: dq,  the result after differentiating dp w.r.t var; 
# 
#=========================================================== 

DF:=proc(dp, var) 
   local dq, dp1, dp2; 
     options remember,system;    

#Step1 [Recursive base] 

      if op(0,dp) <> `+` then 
         dq:=DF_prod(dp,var); 
         RETURN(dq); 
      fi; 
    
#Step2 [Recursion] 

     dp1:=op(1,dp); 
     dp2:=normal(dp-dp1); 
     dq:=normal(DF_prod(dp1,var)+DF(dp2,var));             
          
#Step3 [Prepare for return] 

     RETURN(dq); 
      
end: 

#------------------------------------------ 
# subroutines 
#------------------------------------------ 

#========================================================= 
# 
# 
#                der <- DF_indet(indet, var) 
#       
#           (differentiate a differential indeterminante) 
# 
#   Input: indet, a differential indeterminate which is a 
#                 derivative of symobls in Diff_Indets w.r.t 
#                 some variables in Diff_Var; 
# 
#          var,   the variable w.r.t which to be differeniated 
# 
#   Output: der,  the derivative of indet w.r.t var 
# 
#=========================================================== 

DF_indet:=proc(indet, var) 
         local der, p, i, index, j,Diff_Var,Diff_Indets,Lvar; 
         global depend_relation; 
     options remember,system;          

#STEP0 [Initiation diffs]; 
          
         Diff_Var:=depend_relation[2]; 
         Diff_Indets:=depend_relation[1]; 
         Lvar:=nops(Diff_Var); 
          
#STEP1 [Special cases] 
          
         if  not member(var, Diff_Var, 'p') then 
             der := diff(indet, var); 
             RETURN(der); 
         fi; 
          
         if not member(indet, Diff_Indets) and 
            not member(op(0,indet), Diff_Indets) then 
            der := diff(indet, var); 
            RETURN(der); 
         fi; 
          
#STEP2 [Compute] 

         index:=NULL; #Initialization 
          
         if member(indet, Diff_Indets) then   
          
         # build a derivative from a diff indet       
          
            for i from 1 to Lvar do 
                if i = p then 
                   index:=index,1; 
                else 
                   index:=index,0; 
                fi; 
            od; 
            der:=indet[index]; 
         else   
          
         # modify a derivative from a given one 
          
           for i from 1 to Lvar do 
               j:=op(i,indet); 
               if i = p then 
                  index:=index,(j+1); 
               else 
                  index:=index,j; 
               fi; 
           od; 
           der:=op(0,indet)[index] 
         fi; 
          
#STEP3 [Prepare for return] 
          
          RETURN(der); 
          
end: 

#======================================================== 
# 
#           dq <- DF_power(dp, var) 
# 
#           (differentiating a power of a drivative of a diff indet) 
# 
#           input:  dp, a power of a diff indet 
#                   var, the variable w.r.t which to be differentiated 
# 
#           output: dq, the result after differentiating dp w.r.t var 
# 
#========================================================= 

DF_power:=proc(dp, var) 
         local dq, indet, expon; 
              options remember,system; 
          
#Step1 [Special cases] 

      # constant case 
      
      if dp = 1 then 
         dq := 0; 
         RETURN(dq); 
      fi; 
      
      # indeterminate case 
      
      if op(0,dp) <> `^` then 
         dq := DF_indet(dp, var); 
         RETURN(dq); 
      fi; 
                
#Step2 [Compute] 

      indet:=op(1,dp); 
      expon:=op(2,dp); 
      dq:=expon*indet^(expon-1)*DF_indet(indet,var); 
    
#Step3 [Prepare for return] 

      RETURN(dq); 

end: 

#========================================================= 
# 
#        dq <- DF_prod(dp, var) 
# 
#        (differentiating a product of derivatives of some diff indets)     
# 
#        input: dp, a product of derivatives of some diff indets   
#               var,  the variable w.r.t which to be differentiated 
# 
#        output:  dq, the result after differentiating dp w.r.t var; 
# 
#========================================================== 

DF_prod:=proc(dp, var) 
        local dq, dp1, dp2; 
             options remember,system; 

#Step1 [Recursive base] 

      if op(0,dp) <> `*` then 
         dq := DF_power(dp, var); 
         RETURN(dq); 
     fi; 

#Step2 [Recursion] 

     dp1:=op(1,dp); 
     dp2:=normal(dp/dp1); 
      
     dq:=normal(DF_power(dp1,var)*dp2+dp1*DF_prod(dp2,var)); 
      
#Step3 [Prepare for return] 

     RETURN(dq); 
      
end: 
                                        
#=============================================================== 
# 
#            ord  <- torder(der) 
# 
#           computing the order of a derivative of a diff indet 
# 
#           input: der, a derivative of a diff indet 
# 
#           output: ord, ord is the order of der. 
# 
#================================================================= 

torder:=proc(der) 
local ord,i,Diff_Var,Lvar; 
global depend_relation; 
     options remember,system; 

#STEP1 [Initialize] 

      if type(der,symbol) then RETURN(0) fi; 

         ord := 0; 
              
#STEP2 [Compute] 

           
              for i from 1 to nops(der) do 
                 ord := ord + op(i, der); 
              od; 
 ord; 
end: 
                    

#================================================================== 
#     
#            lead <- Leader(dp, ranking,ord, _cls, _ord) 
#     
#            Input: dp, a  nonzero diff poly; 
#                   ranking, a specified ranking; 
# 
#            Output: the Leader w.r.t. ranking; 
# 
#            Option: _cls, the class of the lead; 
#                   _ord,, the order of the lead; 
# 
#=================================================================== 
Mainvar:=proc(p,ord) 
     options remember,system; 
 Leader(p,ord); 
 end: