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
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
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);
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
434. if not(whattype(nogood)=`+`) then
435.   nogood:=[nogood];
436. else
437.  nogood:=convert(nogood,`list`);
438. fi;
441.
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);