TEXT   15

int

Guest on 16th May 2022 01:50:01 AM

  1. ###############################################################################
  2. #
  3. #  The procedure varder calculates the variational derivative of
  4. #  a functional L with respect to u.
  5. #
  6. ###############################################################################
  7.  
  8. varder:=proc(L,u,x)
  9. option `Authors: B. Deconinck`;
  10. local N,k,lag,v,diffvec,m;
  11. N:=hdifforder(L,u,x);
  12. # N is the order of the highest derivative of u w.r.t. x occuring in L
  13. lag:=L;
  14. for k from N by -1 to 1 do
  15.  lag:=subs(diff(u,[seq(x,j=1..k)])=v[k],lag)
  16. od;
  17. lag:=subs(u=v[0],lag);
  18. diffvec:=[seq(diff(lag,v[k]),k=0..N)];
  19. diffvec:=map2(subs,v[0]=u,diffvec);
  20. diffvec:=map2(subs,[seq(v[k]=diff(u,[seq(x,m=1..k)]),k=1..N)],diffvec);
  21. # The last result serves as the output of the procedure
  22. (expand(diffvec[1]+add((-1)^k*diff(diffvec[k+1],[seq(x,m=1..k)]),k=1..N)))
  23. end:
  24.  
  25. ###############################################################################
  26. #
  27. # The procedure hdifforder finds the order of the highest
  28. # derivative of u w.r.t. x occuring in L.
  29. #
  30. ###############################################################################
  31.  
  32. hdifforder:=proc(L,u,x)
  33. option `Authors: B. Deconinck`;
  34. local k,lijst,flag;
  35. k:=0; flag:=1; lijst:=[];
  36. while flag=1 do
  37.  if has(L,diff(u,lijst)) then
  38.   lijst:=[op(lijst),x];
  39.   k:=k+1
  40.  else flag:=0
  41.  fi
  42. od;
  43. # the last result serves as the output of the procedure
  44. if k=0 then 0
  45. else k-1
  46. fi;
  47. end:
  48.  
  49. ###############################################################################
  50. #
  51. # The procedure varder calculates the higher-order
  52. # variational derivative of a functional L with respect to u.
  53. # This procedure cannot be used to compute the variational
  54. # derivative. We require n to be at least 1.
  55. #
  56. ###############################################################################
  57.  
  58. vardern:=proc(L,u,n,x)
  59. option `Authors: B. Deconinck`;
  60. local N,k,lag,v,diffvec,m;
  61. N:=hdifforder(L,u,x);
  62. if N<n then return(0)
  63. fi;
  64. # N is the order of the highest derivative of u w.r.t. x occuring in L
  65. lag:=L;
  66. for k from N by -1 to n do
  67.  lag:=subs(diff(u,[seq(x,j=1..k)])=v[k],lag)
  68. od;
  69. diffvec:=[seq(diff(lag,v[k]),k=n..N)];
  70. diffvec:=map2(subs,[seq(v[k]=diff(u,[seq(x,m=1..k)]),k=n..N)],diffvec);
  71. # The last result serves as the output of the procedure
  72. expand(diffvec[1]+add((-1)^(k-n)*binomial(k,n)*diff(diffvec[k-n+1],[seq(x,m=n..k-1)]),k=n+1..N))
  73. end:
  74.  
  75. ################################################################################################
  76. #
  77. # This is a modified homotopy integral, a list of terms is returned
  78. # Here u is a vector of unknown functions
  79. # that depend on an independent variable x.
  80. #
  81. ################################################################################################
  82.  
  83. vechomotopyintegral:=proc(LL,u,x)
  84. option `Authors: B. Deconinck, M. Nivala`;
  85. local sh,result,newans,newflux,aa,limset,ans,i,r,m,lambda,lambda0,flux,scaledflux,summedflux,N,n,L,h,a,s;
  86. L:=expand(eval(LL));
  87. n:=nops(u); # number of unknown functions
  88. # N is the vector of differential orders of u in L
  89. N:=[seq(hdifforder(L,u[r],x),r=1..n)];
  90. # flux is the integrand before scaling of the arguments
  91. # to obtain the maximum number of terms, we take this as a list and a sum
  92. flux:=expand([seq(op([u[r]*vardern(L,u[r],1,x),seq(diff(u[r]*vardern(L,u[r],i+1,x),[seq(x,m=1..i)]),i=1..N[r]-1)]),r=1..n)]);
  93. scaledflux:=expand(subs({seq(u[r]=lambda*u[r],r=1..n)},flux));
  94. sh:=normal(convert(scaledflux,`+`));
  95. ans:=expand({seq(int(expand(scaledflux[i]/lambda),lambda),i=1..nops(scaledflux)),int(expand(sh/lambda),lambda)});
  96. #remove expressions that were not integrated
  97. for s in ans do
  98.  if has(s,int) then
  99.   ans:=ans minus {s}
  100.  fi;
  101. od;
  102. #get as many terms as possible
  103. newans:=[NULL];
  104. for i from 1 to nops(ans) do
  105.  if whattype(ans[i])=`+` then
  106.   newans:=[op(newans),op(ans[i])];
  107.  else
  108.   newans:=[op(newans),ans[i]];
  109.  fi;
  110. od;
  111. #normalize the solution and remove constants
  112. result:={NULL};
  113. for i from 1 to nops(newans) do
  114.  if has(newans[i],x) then
  115.   result:={op(result),SPLIT(newans[i],x)[1]};
  116.  fi;
  117. od;
  118. [op(subs(lambda=1,result))];
  119. end:
  120.  
  121. ###############################################################################
  122. #
  123. # whatfun picks out undefined functions in an expression
  124. #
  125. ###############################################################################
  126.  
  127. whatfun:=proc(E)
  128. option `Authors: M. Nivala`;
  129. local n,s,sid;
  130. s:={NULL};
  131. sid:=[op(indets(E,`function`))];
  132. for n from 1 to nops(sid) do
  133.  if not(has(sid[n],`diff`)) and not(type(op(0,sid[n]),`procedure`))
  134.   then s:={op(s),sid[n]};
  135.  fi;
  136. od;
  137. [op(s)];
  138. end:
  139.  
  140. ##############################################################################
  141. #
  142. # orderset orders a set from large to small
  143. #
  144. ##############################################################################
  145.  
  146. orderset:=proc(s)
  147. option `Authors: B. Deconinck, M. Nivala`;
  148.     local ss,leftover,k;
  149.     ss:=NULL;
  150.     leftover:=s;
  151.     for k to nops(s) do
  152.         ss:=ss,max(op(leftover));
  153.         leftover:=leftover minus {max(op(leftover))}
  154.     od;
  155.     [ss]
  156. end:
  157.  
  158. #################################################################################
  159. #
  160. # funlexorderconverts a list of functional expressions to strings, orders them
  161. # according to lexorder, and returns the funtional expressions
  162. #
  163. #################################################################################
  164.  
  165. funlexorder:=proc(LIST)
  166. option `Authors: M. Nivala`;
  167. local Newlist,i,subrules,s;
  168. subrules:=[seq(convert(LIST[i],`string`)=LIST[i],i=1..nops(LIST))];
  169. Newlist:=[seq(op(1,subrules[i]),i=1..nops(LIST))];
  170. s:=sort(Newlist,`lexorder`);
  171. subs(subrules,s);
  172. end:
  173.  
  174. #######################################################################################
  175. #
  176. # orders a set according to order of derivatives, from highest to lowest,
  177. # grouping those of the same order
  178. #
  179. #######################################################################################
  180.  
  181. dorderlist:=proc(L,u,x)
  182. option `Authors: M. Nivala`;
  183. local finlist,j,diffset,dlist,i,MAX,difflist,ord;
  184. difflist:=[seq(hdifforder(L[i],u,x),i=1..nops(L))];
  185. diffset:=orderset({op(difflist)});
  186. finlist:=[NULL];
  187. for i in diffset do
  188.  dlist:=[NULL];
  189.  for j from 1 to nops(L) do
  190.   if i=difflist[j] then
  191.    dlist:=[op(dlist),L[j]];
  192.   fi;
  193.  od;
  194. finlist:=[op(finlist),dlist];
  195. od;
  196. finlist;
  197. end:
  198.  
  199. ###################################################################################
  200. #
  201. # torderset orders a set according to number of products in a term
  202. #
  203. ###################################################################################
  204.  
  205. torderset:=proc(E)
  206. option `Authors: M. Nivala`;
  207. local w,fin,n,s;
  208. s:={seq(nops(op(1,E[n]))+(n/(nops(E)+1)),n=1..nops(E))};
  209. s:=orderset(s);
  210. w:={seq(nops(op(1,E[n]))+(n/(nops(E)+1))=E[n],n=1..nops(E))};
  211. fin:=subs(w,s);
  212. end:
  213.  
  214. ###################################################################################
  215. #
  216. # SPLIT splits the product T into an x-dependent part and a constant part
  217. #
  218. ###################################################################################
  219.  
  220. SPLIT:=proc(T,x)
  221. option `Authors: M. Nivala`;
  222. local xpart,coeffpart,s,opT,xset,coeffset;
  223. xset:={NULL};
  224. coeffset:={NULL};
  225. if whattype(T)=`*` then
  226.  opT:=[op(T)];
  227. else
  228.  opT:=[T];
  229. fi;
  230. for s in opT do
  231.  if has(s,x) then
  232.   xset:={op(xset),s};
  233.  else
  234.   coeffset:={op(coeffset),s};
  235.  fi;
  236. od;
  237. xpart:=convert(xset,`*`);
  238. coeffpart:=convert(coeffset,`*`);
  239. return([xpart,coeffpart]);
  240. end:
  241.  
  242. #################################################################################
  243. #
  244. # SUBS substitutes ssubrules into EE
  245. #
  246. #################################################################################
  247.  
  248. SUBS:=proc(ssubrules,EE,x)
  249. option `Authors: M. Nivala`;
  250. local E,Elist,subrules;
  251. subrules:=torderset(ssubrules);
  252. E:=expand(EE);
  253. if not(whattype(E)=`+`) then
  254.  Elist:=[E];
  255. else
  256.  Elist:=convert(E,`list`);
  257. fi;
  258. Elist:=[seq(SPLIT(Elist[i],x),i=1..nops(Elist))];
  259. Elist:=subs(subrules,Elist);
  260. E:=add(convert(Elist[i],`*`),i=1..nops(Elist));
  261. end:
  262.  
  263. ##################################################################################
  264. #
  265. # GROUP goups terms with the same x-dependence in E
  266. #
  267. ##################################################################################
  268.  
  269. GROUP:=proc(E,x)
  270. option `Authors: M. Nivala`;
  271. local EL,subrules,EE,Esubs,newsubrules,i,b;
  272. if not(whattype(E)=`+`) then
  273.   EL:=[E];
  274. else
  275.  EL:=convert(E,`list`);
  276. fi;
  277. subrules:=[op({seq(SPLIT(EL[i],x)[1]=b[i],i=1..nops(EL))})];
  278. subrules:=torderset(subrules);
  279. Esubs:=SUBS(subrules,E,x);
  280. Esubs:=collect(Esubs,{seq(b[i],i=1..nops(subrules))});
  281. newsubrules:=[seq(op(2,subrules[i])=op(1,subrules[i]),i=1..nops(subrules))];
  282. EE:=subs(newsubrules,Esubs);
  283. end:
  284.  
  285. ################################################################################################
  286. #
  287. # RRED performs gaussian elimination with a minimal number of row switching operations
  288. # it row-reduces A while preserving the order of the rows
  289. # returns the first maxsstep-1 linearly independent rows
  290. #
  291. ################################################################################################
  292.  
  293. RRED:=proc(A,sstep,B,maxsstep)
  294. option `Authors: M. Nivala`;
  295. local step,i,j,newA,flag1,rowdims,newB;
  296. rowdims:=linalg[rowdim](A);
  297. if sstep=maxsstep then
  298.  return(B);
  299. fi;
  300. step:=sstep;
  301. newA:=A;
  302. flag1:=0;
  303. for i from 1 to rowdims do
  304.   if newA[i,step]<>0 and flag1=0 then
  305.    flag1:=i;
  306.   fi;
  307. od;
  308. if flag1<>0 then
  309.  newA:=linalg[pivot](newA,flag1,step);
  310.  if step=1 then
  311.   newB:=linalg[matrix]([linalg[row](newA,flag1)]);
  312.  else
  313.   newB:=linalg[stackmatrix](B,linalg[row](newA,flag1));
  314.  fi;
  315.  if rowdims=1 then
  316.   return(newB);
  317.  else
  318.   newA:=linalg[delrows](newA,flag1..flag1);
  319.   RRED(newA,step+1,newB,maxsstep);
  320.  fi;
  321. else
  322.  RRED(newA,step+1,newB,maxsstep);
  323. fi;
  324. end:
  325.  
  326. #################################################################################################################################
  327. #
  328. # INT is used to compute integrals of expressions that depend on unknown functions that depend on an independent variable x.
  329. # All terms that are total derivatives are integrated.
  330. # The order of the highest derivative in the remaining integral is minimized.
  331. #
  332. #################################################################################################################################
  333.  
  334. INT:=proc(EE,x)
  335. option `Authors: M. Nivala`;
  336. local dcheck,limset,eqnset2,RE,param, eqnset,v,slnset,novan,Flist,intnogood,nogood,subDFL,intxset,sln,Esub,Fsub,DFLsub,termlist,FL,xset,E,u,goo,A,Ar,Ag,bob,syst,i,j,m,DFL,s,n,fresult,result,a,E1,DF1,R,subrules,b,EL,F1,F,DF,k;
  337.  
  338. E:=expand(eval(EE));
  339.  
  340. #determine the unknown functions
  341. u:=whatfun(E);
  342. if u=[NULL] then return(int(E,x))
  343. fi;
  344.  
  345. #check if E is exact
  346. #dcheck:=[seq(varder(E,u[i],x),i=1..nops(u))];
  347. #if dcheck=[seq(0,i=1..nops(u))] then
  348.  #return(expand(int(E,x)));
  349. #fi;
  350.  
  351. #convert E to a list
  352. if not(whattype(E)=`+`) then
  353.  EL:=[E];
  354. else
  355.  EL:=convert(E,`list`);
  356. fi;
  357.  
  358. #remove terms which have no functional or x dependence to be integrated directly
  359. xset:={NULL};
  360. for k in EL do
  361. if not(has(k,u)) or not(has(k,x)) then
  362.  xset:={op(xset),k};
  363. fi;
  364. od;
  365. xset:=expand(convert(xset,`+`));
  366. intxset:=int(xset,x);
  367. E:=expand(eval(E-xset));
  368. if not(whattype(E)=`+`) then
  369.  EL:=[E];
  370. else
  371.  EL:=convert(E,`list`);
  372. fi;
  373.  
  374. #FL is a list of terms from the homotopy integral of E
  375. FL:=vechomotopyintegral(E,u,x); #print(FFFFFFF,FL);
  376. if FL=[NULL] then
  377.  return(int(E,x)+intxset);
  378. fi;
  379.  
  380. #DFL is the derivative of FL with respect to x
  381. DFL:=[seq(expand(diff(FL[n],x)),n=1..nops(FL))];
  382.  
  383. #form substitution rules
  384. #subDFL is for substitution rules
  385. subDFL:={NULL};
  386. for k in DFL do
  387.  if not(whattype(k)=`+`) then
  388.    subDFL:={op(subDFL),k};
  389.  else
  390.   subDFL:={op(subDFL),op(k)};
  391.  fi;
  392. od;
  393. termlist:=[op({op(EL),op(FL),op(subDFL)})];
  394. #remove constant coefficients
  395. termlist:=[op({seq(SPLIT(termlist[i],x)[1],i=1..nops(termlist))})];
  396. #orders termlist according to the order of the highest derivative of each element, from highest to lowest, groouping those of the
  397. #same order
  398. termlist:=dorderlist(termlist,u,x);#print(dord,termlist);
  399. #for consistency, implement an ordering of funtions and their derivatives according to lexorder for those terms which are of the
  400. #same differential order
  401. termlist:=[seq(op(funlexorder(termlist[i])),i=1..nops(termlist))];#print(orderedsub,termlist);
  402.  
  403. #perform the substitutions
  404. subrules:=(seq(termlist[n]=b[n],n=1..nops(termlist)));
  405. #order subrules according to number of products in each term
  406. subrules:=torderset({subrules});#print(SUBRULES,subrules);
  407. Esub:=SUBS(subrules,E,x);
  408. DFLsub:=[seq(SUBS(subrules,DFL[i],x),i=1..nops(DFL))];
  409.  
  410. #R is what remains after integrating E
  411. #R is optimized, removing higest order derivaites first
  412. R:=expand(Esub-sum(a[m]*DFLsub[m],m=1..nops(DFL)));#print(E-sum(a[m]*DFL[m],m=1..nops(DFL)));
  413. R:=collect(R,{seq(b[m],m=1..nops(termlist))});#print(RRR,R);
  414. #form the linear system to be solved
  415. syst:={NULL};
  416. for m from 1 to nops(termlist) do
  417.  if has(coeff(R,b[m]),[seq(a[n],n=1..nops(DFL))]) then
  418.   syst:=[op(syst),coeff(R,b[m])=0];
  419.  fi;
  420. od;
  421.  
  422. #A is the overdetermined system for the coefficients, a
  423. #Ar is the first nops(a) linearly independent rows of A
  424. A:=linalg[genmatrix](syst,[seq(a[j],j=1..nops(DFLsub))],flag);#print(AAAA,A);
  425. Ar:=RRED(A,1,A,nops(DFLsub)+1);#print(ArAr,Ar);
  426.  
  427. slnset:=linalg[backsub](Ar,`false`,`v`);#print(slnset);
  428. sln:=seq(a[m]=slnset[m],m=1..nops(DFLsub));#print(solution,sln);
  429. sln:=eval(sln,v=0);#print(solution,sln);
  430. assign(sln);
  431.  
  432. #all terms which were not integrated using homotopy, attempt to integrate each independently without homotopy
  433. nogood:=expand(E-add(a[i]*DFL[i],i=1..nops(DFL)));#print(nogood);
  434. if not(whattype(nogood)=`+`) then
  435.   nogood:=[nogood];
  436. else
  437.  nogood:=convert(nogood,`list`);
  438. fi;
  439. intnogood:=expand(add(int(nogood[i],x),i=1..nops(nogood)));#print(intnogood);
  440. #print(expand(E-add(a[i]*DFL[i],i=1..nops(DFL))+xset+novan+RE),"Was Integrated By Maple");
  441.  
  442. fresult:=add(a[j]*FL[j],j=1..nops(FL))+combine(intnogood)+intxset;#print(fdfdfdfd,fresult);
  443. unassign(a);
  444.  
  445. #check if correct
  446. if simplify(expand(EE-diff(fresult,x)))<>0 then
  447.  print("No integration done.");
  448.  return(int(E,x)+intxset)
  449. else
  450.  return(fresult);
  451. fi;
  452. end:
  453.  
  454. ########################################################################################################
  455. #
  456. # ouputs the expression G modulo total derivatives
  457. #
  458. ########################################################################################################
  459.  
  460. DerivativeMod:=proc(G,x)
  461. local IG,n,intpart,Glist;
  462. IG:=INT(G,x);
  463. if not(whattype(IG)=`+`) then
  464.  Glist:=[IG];
  465. else
  466.  Glist:=convert(IG,`list`);
  467. fi;
  468. intpart:=0;
  469. for n from 1 to nops(Glist) do
  470.  if has(Glist[n],`int`) then intpart:=diff(Glist[n],x);
  471.  fi;
  472. od;
  473. intpart;
  474. end:
  475.  
  476. #########################################
  477. #
  478. # checks if INT is correct
  479. #
  480. #########################################
  481.  
  482. INTcheck:=proc(E,F,x)
  483. normal(expand(E-diff(F,x)));
  484. end:
  485.  
  486. ##############################################
  487. #
  488. #
  489. #
  490. ##############################################
  491.  
  492. chi:=proc(n,x)
  493. option remember;
  494. local result;
  495. if n=0 then result:=-u(x);
  496.  else if n=1 then result:=diff(u(x),x);
  497.          else result:=-diff(chi(n-1,x),x)-add(chi(n-k-2,x)*chi(k,x),k=0..n-2);
  498.       fi;
  499. fi;
  500. result;
  501. end:

Raw Paste


Login or Register to edit or fork this paste. It's free.