Re: For Loop and NDSolve

*To*: mathgroup at smc.vnet.net*Subject*: [mg85209] Re: For Loop and NDSolve*From*: Szabolcs <szhorvat at gmail.com>*Date*: Sat, 2 Feb 2008 03:25:41 -0500 (EST)*References*: <fnuh4s$9ni$1@smc.vnet.net>

On Feb 1, 8:18 am, PV <prasanna.p... at gmail.com> wrote: > I am trying to solve a time dependent ordinary differential equation > of matrices. Here my coefficient matrices in the DE change at every > instant of time, depending on the solution of the DE in the previous > instant( feed back type problem). To implement this I use NDSolve with > a For loop in time. This works fine for certain size of my matrices > but has got me into trouble for other sizes. I run out of memory > pretty quickly. I have 1 GB Ram with an AMD X2 duo core 1.7GHZ > processor. > Here are the routines that cause problems:: > > N0 = 1921(large number) > > For[j=1,j\[LessEqual]100, > t = j*delt; > solnmatrix = > SparseArray[{{i_,i_}\[Rule](-V0[a]/2 + ((i-(4*N0+1))*k0-2*t)^2), > {m_,n_}/; > Abs[m-n]\[Equal]N0\[Rule]-V0[a]/4},{24*N0+1,24*N0+1}] ; > Print["pbegin ",MemoryInUse[]]; > soln = NDSolve[{I*c3*u'[p]== solnmatrix.u[p],u[0]\[Equal]wavefn}, > u,{p,0,delt}]; > Print["p1 ",MemoryInUse[]]; > wavefn = First[u[delt]/.soln]; > Clear[solnmatrix]; > Clear[soln]; > Print["p2 ",MemoryInUse[]]; > C1 = ( Tr[Abs[wavefn]^2/2]+0.25* > Sum[Conjugate[wavefn[[n]]]*wavefn[[n+N0]],{n,1,23*N0+1}] + > 0.25*Sum[Conjugate[wavefn[[n]]]*wavefn[[n-N0]], > {n,N0+1,24*N0+1}])* > L; > b= -\[Eta]/(I*\[Gamma]* > C1 - \[Kappa]) + ((\[Eta] + (I*\[Gamma]*C1 - \[Kappa])* > a )/(I*\[Gamma]*C1 - \[Kappa]))* > E^(delt*Tb*(I*\[Gamma]*C1 - \[Kappa])); > a = b; > AppendTo[fieldlist3,a]; > AppendTo[normlist,L*Tr[Abs[wavefn]^2]]; > AppendTo[couplinglist,C1]; > j++; > Print["pend ",MemoryInUse[]]; > ] > > One way I tried to address this issue is by defining the operation > that I iterate inside my For loop as a function, and then look at how > much memory each evaluation of that function consumes. > The function is:: > timevol[a_,j_]:= > Block[{t=j*delt,wave=wavefn,soln,solnmatrix,C1}, > solnmatrix = > SparseArray[{{i_,i_}\[Rule](-V0[a]/2 + ((i-(4*N0+1))*k0-2*t)^2), > {m_,n_}/; > Abs[m-n]\[Equal]N0\[Rule]-V0[a]/4},{24*N0+1,24*N0+1}] ; > soln = > NDSolve[{I*c3*u'[p]== solnmatrix.u[p],u[0]\[Equal]wave},u,{p, > 0,delt}]; > Share[]; > wave = First[u[delt]/.soln]; > wavefn = wave; > C1 = ( > Tr[Abs[wavefn]^2/2]+0.25* > Sum[Conjugate[wavefn[[n]]]*wavefn[[n+N0]],{n,1,23*N0+1}] > + > 0.25*Sum[Conjugate[wavefn[[n]]]*wavefn[[n-N0]], > {n,N0+1,24*N0+1}])* > L; > b= -\[Eta]/(I*\[Gamma]* > C1 - \[Kappa]) + ((\[Eta] + (I*\[Gamma]*C1 - \ > [Kappa])* > a )/(I*\[Gamma]*C1 - \[Kappa]))* > E^(delt*Tb*(I*\[Gamma]*C1 - \[Kappa])); > Clear[solnmatrix]; > Clear[soln]; > Unprotect[In,Out]; > Clear[In,Out]; > b > ] > > The memory consumption of this function is also very high and it > changes every time I evaluate it, espescially when j value is large it > starts using more memory. Hi, While I did not try to untangle the calculations (the code does not run on its own, as you sent it), here are a few tips for making the code more manageable, and possibly more efficient: - Never ever use For in Mathematica! Especially for building up an array of the results. Instead of this: arr = {} For[i = 1, i <= 10, i++, AppendTo[arr, i^2] ] Use this: Table[i^2, {i, 1, 10}] - Do not use Clear[] inside functions/loops (such as Table[]). For local/temporary variables, use a Module[] (it is better to use a Module[] instead of a Block[] until you understand the difference between them, though it won't matter in this particular case). This will ensure that no temporary values will escape from the calculation to pollute the Global context and eat up memory. When possible, consider using With[]. - Try to avoid explicit indexing and work directly with vectors/array whenever possible. Use Total@Conjugate[ wavefn[[1 ;; 23 N0 + 1]] wavefn[[N0 ;; 24 N0 + 1]] ] instead of Sum[Conjugate[wavefn[[n]]]*wavefn[[n + N0]], {n, 1, 23*N0 + 1}] To be honest, I don't think that these changes will help a lot with memory usage (getting rid of the AppendTo[]s might help a little bit). Some operations simply require a lot of memory. But at least they will help keep things clean and manageable and make it easier to figure out where exactly the memory is used up. In/Out are not being filled up during a calculation, so don't clear them, especially not inside a function. If you are really worried about them, evaluate $HistoryLength = 0 right after the kernel starts up.