using structures and iterators
- To: mathgroup at yoda.physics.unc.edu
- Subject: using structures and iterators
- From: rclark at lpl.arizona.edu (Richard Clark x4971)
- Date: Thu, 7 Apr 94 07:47:33 MST
I am trying to find the propper use of Mma iterators and object structures
to simplify a proceedure. I have two equivalent sets of equations and I
need to find the relation between the sets of variables in which they are
expressed to a desired order of aproximation.
The variables being solved for are related to the extent that var1 is
first order, var2 is 2nd order ... so that in the lower order partial
solution the higher order varN (of either set) do not appear.
x1,y1 ~m largest term in eqn1 is ~m
x2,y2 ~m^2 largest term in eqn2 is ~m^2
x3,y3 ~m^3 largest term in eqn3 is ~m^3
and m << 1
I find the solutions to successively higher order:
solve eqn1 for x1 to order m
substitute this in eqn2 for and solve x2 to order m^2
substitute this in eqn1 for and solve x1 again, now to order m^2
substitute this in eqn3 for and solve x3 to order m^3
substitute this in eqn2 for and solve x2 again, now to order m^3
substitute this in eqn1 for and solve x1 again, now to order m^3
At present I am working only to 3rd order (which has already had a partial
solution published). The following code works by performing each of the
above steps explicitly. I want to extend this solution to higher order
accuracy of solution for {x1,x2,x3,x4...}. Also this same type of
successive aproximation solution will be used in several different steps
in obtaining the complete solution to the response functions I am working
on.
Any suggestions would be apreciated.
Richard Clark
rclark at lpl.arizona.edu
------------------------
Now for the listings. (These variable names follow some traditional
notation for the problem.)
(e, Jj[1]) 1st order (Jj[1] is sometimes replaced with another 1st order
(k, Jj[2]) 2nd order quantity 'ep' for reasons that were convienent
(h, Jj[3]) 3rd order half a century ago)
This first one has the advantage that it *does* work.
I want to generalize it by using Mma's iterators instead of an editor!
--------------------------------------------------------
(* hkslv: solve (ep,Jn) <--> (e,k,h) *)
(* get the 2 equivalent sets of equations *)
<<rslt/ra.res
<<rslt/rahk.res
(* define the order of the terms *)
qhrul={q->aa q,ep->aa ep,e->aa e,k->aa^2 k,h->aa^3 h,Jj[n_]->aa^n Jj[n]}
ehkrul={ee[n_]->aa^n ee[n], kk[n_]->aa^(n+1) kk[n], hh[n_]->aa^(n+2) hh[n]}
(* after substitution express e,k,h as sums of various orders *)
ehksub={e->ees, k->kks, h->hhs}
ees=Sum[ee[n],{n,1,maxord}]
kks=Sum[kk[n],{n,1,maxord-1}]
hhs=Sum[hh[n],{n,1,maxord-2}]
(* make copies of the 2 systems containing terms up to a desired order *)
Subser[ord_] := Block[{i},
Do[sshktmp[i]=Normal[Series[sshk[i]/.ehksub/.qhrul/.ehkrul,\
{aa,0,ord}]/.aa->1],{i,1,maxord}];
Do[sstmp[i]=Normal[Series[ss[i]/.qhrul,\
{aa,0,ord}]/.aa->1],{i,1,maxord}];
sstm=Table[sstmp[i]==sshktmp[i],{i,1,maxord}];
]
(* 1st order *)
Subser[1]
e1=Solve[sstm[[1]],ee[1]]
ee[1]=ee[1]/.e1[[1]]
(* 2nd order *)
Subser[2]
k1=Solve[sstm[[2]],kk[1]]
kk[1]=kk[1]/.k1[[1]]
e2=Solve[sstm[[1]],ee[2]]
ee[2]=ee[2]/.e2[[1]]
(* 3rd order *)
Subser[3]
h1=Solve[sstm[[3]],hh[1]]
hh[1]=hh[1]/.h1[[1]]
k2=Solve[sstm[[2]],kk[2]]
kk[2]=kk[2]/.k2[[1]]
e3=Solve[sstm[[1]],ee[3]]
ee[3]=ee[3]/.e3[[1]]
ees=Expand[ees]
kks=Expand[kks]
hhs=Expand[hhs]
Save["hkslv.res", ees, kks, hhs]
--------------------------------------------------------
This is what I've come up with so far. It looks like it should be
close to what I want but doesn't work.
(* hkslv: solve (ep,Jn) <--> (e,k,h) (21-30) *)
<<rslt/ra.res
<<rslt/rahk.res
maxord=3
(* define the order of terms *)
qhrul={q->aa q,ep->aa ep,e->aa e,k->aa^2 k,h->aa^3 h,Jj[n_]->aa^n Jj[n]}
(* make a list of the variables to solve for *)
ehk={e, k, h}
(* called as sln=Solser[3, sshk, ss, ehk] *)
(* set1 is the system expressed in varset *)
Solser[ord_,set1_,set2_,varset_]:=Block[{},
vsub=Table[varset[[i]]->v[i],{i,1,maxord}];
vrul={vn[n_,m_]->aa^m vn[n,m]};
(* ideally the solution will end up in v *)
(* I count the order of the vn differently from [ehk]s in version 1
so that the second subscript is actually the order of the term *)
Do[v[i]=Sum[vn[i,j],{j,i,maxord}],{i,1,maxord}]
Do[Block[{},
(* these 3 statements are the equivalent of Subser in version 1 *)
Do[s1s[i]=Normal[Series[set1[i]/.vsub/.qhrul/.vrul,\
{aa,0,ord}]/.aa->1],{i,1,iord}];
Do[s2s[i]=Normal[Series[set2[i]/.qhrul,\
{aa,0,ord}]/.aa->1],{i,1,iord}];
stm=Table[s1s[i]==s2s[i],{i,1,iord}];
Do[Block[{},
slv=Solve[stm[[iord-jord+1]],vn[iord-jord+1,iord]];
vn[iord-jord+1,iord]=vn[iord-jord+1,iord]/.slv[[1]];
],{jord,1,iord}];
],{iord,1,ord}]
]
--------------------------------------------------------
Here is ./rslt/ra.res
maxord = 3
ss/: ss[0] = (-4*ep^2)/45 - (512*ep^3)/2835 + (4*ep^2*q)/35
ss/: ss[1] = (-2*ep)/3 - (44*ep^2)/63 - (1286*ep^3)/945 + (3*ep*q)/7 +
(25*ep^2*q)/21 - (4*ep*q^2)/21 - (9*ep*Jj[2])/14 - (25*q*Jj[2])/168
ss/: ss[2] = (-16*ep^2)/35 - (1376*ep^3)/1155 + (4*ep*q)/7 +
(236*ep^2*q)/231 - (18*ep*q^2)/77 - Jj[2] - (548*ep*Jj[2])/231 -
(18*q*Jj[2])/77
ss/: ss[3] = (-16*ep^3)/33 + (160*ep^2*q)/231 - (100*ep*q^2)/231 -
(20*ep*Jj[2])/11 + (25*q*Jj[2])/33 - Jj[3]
--------------------------------------------------------
and ./rslt/rahk.res
maxord = 3
sshk/: sshk[0] = (-4*e^2)/45 - (52*e^3)/567 - (32*e*k)/315
sshk/: sshk[1] =
(-2*e)/3 - (23*e^2)/63 - (4*e^3)/27 + (2*h)/21 - (8*k)/21 - (152*e*k)/315
sshk/: sshk[2] =
(12*e^2)/35 + (4*e^3)/11 + (192*h)/385 + (32*k)/35 + (32*e*k)/105
sshk/: sshk[3] = (-40*e^3)/231 - (80*h)/231