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.