Re: Iterative Type Programming
- To: mathgroup@smc.vnet.net
- Subject: [mg11172] Re: Iterative Type Programming
- From: Allan Hayes <hay@haystack.demon.co.uk>
- Date: Wed, 25 Feb 1998 03:31:44 -0500
- References: <6ctaqj$l4b@smc.vnet.net>
Corey & Alicia Beaverson wrote: > > Below is a program that I wrote as the solution to the temperatures > at different nodes in a "gridwork" (a homework problem). I am still not > as proficient with Mathematica as I would like to be and would like > some opinions on a better way to write an iterative type program. One > of the requirements was to always use the latest values when solving > the equations (Gauss-Seidel). I really feel like I fumbled through this > thing! Also, at the time the only way I could get it to converge was by > making it iterate until "maxiter" was met, but I wanted to have it > iterate until the "tol" was met, any suggestions on how to do this as > well? > > Off[General::spell]; gs := > Tn = Array[Tnew, 4]; > To = Array[Told, 4]; > error = Array[err, 4]; > iter = 0; > ji = 10; > hi = 15; > r = 30; > s = 35; > t = 50; > Tnew[1] = 100; Tnew[2] = 900; Tnew[3] = 1.7; Tnew[4] = 30; Told[1] = 0; > Told[2] = 0; Told[3] = 0; Told[4] = 0; maxiteration = 100; > Tol = 0.001; > While[iter < maxiteration, iter++; If[maxerr > Tol, To = Tn]; > Tn[[1]] = (ji*hi + To[[2]] + To[[3]])/4; > Tn[[2]] = (ji*r + To[[1]] + To[[4]])/4; > Tn[[3]] = (ji*s + To[[1]] + To[[4]])/4; > Tn[[4]] = (ji*t + To[[2]] + To[[3]])/4]; > Do[err[i] = Abs[Told[i] - Tnew[i]], {i, 1, 4}]; maxerr = Max[error]; > On[General::spell]; > Print[TableForm[N[Tn]]] I hope that the following will help. Usually with Mathematica we aim at making a function to do the job: gs[T_,K_, maxiter_,tol_] = procedure. The test inputs will be In[1]:= K = 10 {15.,30.,35.,50.}; T = {100.,900.,1.7,30.}; Use lists and list operations. Localise variables used in the procedure inside Module Calculate error *inside* the While loop. In[2]:= gs[T_,K_, maxiter_,tol_]:= Module[{To=0 T,Tn = T,iter =0,Ton}, While[ iter < maxiter && Max[Abs[To-Tn]]>tol, iter++; Ton = Tn; (*next To*) Tn = (K+ To[[{2,1,1,2}]] + To[[{3,4,4,3}]])/4; To=Ton (*set new To*) ]; Tn (*output value of Tn*) ] In[3]:= gs[T,K,100,.001] Out[3]= {118.751,156.25,168.75,206.251} The extra variable Ton can be avoided by more use of lists (as {a,b} = {1,2}) In[4]:= gs[T_,K_, maxiter_,tol_]:= Module[{To=0 T,Tn = T,iter =0}, While[ iter < maxiter && Max[Abs[To-Tn]]>tol, iter++; {To ,Tn }= {Tn, (K+ To[[{2,1,1,2}]] + To[[{3,4,4,3}]])/4} ]; Tn ] In[5]:= gs[T,K,100,.001] Out[5]= {118.751,156.25,168.75,206.251} FixedPoint version The function FixedPoint is designed to do this sort of thing In[6]:= gs[T_,K_, maxiter_,tol_]:= FixedPoint[ (K+ #[[{2,1,1,2}]] + #[[{3,4,4,3}]])/4&,(*a pure function*) T,maxiter, SameTest -> (Max[Abs[#1-#2]]<=tol&) ] In[7]:= gs[T,K,100,.001] Out[7]= {118.75,156.25,168.75,206.25} It is often useful to display the dynamics by defining helper functions: With In[8]:= update[T_] := (K+ T[[{2,1,1,2}]] + T[[{3,4,4,3}]])/4 error[x_,y_] := Max[Abs[x-y]] We get In[10]:= gs[T_,K_, maxiter_,tol_]:= FixedPoint[update,T,maxiter, SameTest -> (error[#1,#2]<=Tol&)] In[11]:= gs[T,K,100,.001] Out[11]= {118.75,156.25,168.75,206.25} -- Allan Hayes Mathematica Training and Consulting Leicester, UK hay@haystack.demon.co.uk http://www.haystack.demon.co.uk voice: +44 (0)116 271 4198 fax: +44 (0)116 271 8642