For Loop and NDSolve

*To*: mathgroup at smc.vnet.net*Subject*: [mg85192] For Loop and NDSolve*From*: PV <prasanna.phys at gmail.com>*Date*: Fri, 1 Feb 2008 02:16:58 -0500 (EST)

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. It would be great if I could get a solution for this!! Cheers Prasanna