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