Time, Inverse, Simulation simple dynamical system, speed issue

*To*: mathgroup at smc.vnet.net*Subject*: [mg77380] Time, Inverse, Simulation simple dynamical system, speed issue*From*: kristoph <kristophs.post at web.de>*Date*: Thu, 7 Jun 2007 04:03:33 -0400 (EDT)

Attached you find a simulation of 5 dynamical systems each consisting of 1000 Periods. For each period t = 1,...,1000 an inverse of a 2x2m matrix needs to be computed. It takes about 6 seconds to simulate the 5 dynamical systems. I would be grateful for any hint in getting the simulation done much quicker. It seems that calculating the inverse Wel = N[Inverse[ IdentityMatrix[ 2] - ThetaMatrix.LamRho, Method -> \ DivisionFreeRowReduction].ThetaMatrix.divi, 1000] takes most of the time, therefore I tried to solve the problem numerically (see above). This is a lot faster then using just Inverse[...]. I tried LinearSolve, but it is slower then the above. Decreasing the precision can result in a significant error, depending on the parameter constellation. As you might have guessed I have to simulate not only 5 dynamical systems for different parameter constellations. I would also like to increase the order of the matrix that needs to be inverted. Here is the code. Thanks in advance. Kristoph << Statistics`ContinuousDistributions`;(*packages and starting values needed for the simulation*) << Graphics`MultipleListPlot`; Wel = {5, 5}; Lam = {{1/2, 2/3}, {1/2, 1/3}}; Rho = {{1/2, 0}, {0, 648/1000}}; RelDiv = {{2/3, 1/3}, {1/3, 2/3}, {1, 0}}; theta[i_, k_] := Lam[[i, k]]*Rho[[i, i]]*Wel[[i]]/Sum[Lam[[j, k]]*Rho[[j, j]]*Wel[[j]], {j, 1, 2}]; Timing[For[n = 1, n =E2=89=A4 5, n++,(*begin loop for the 5 simulations o= f the dynamical systeme*) Wel = {5, 5};(*starting values*) RelWel = {};(*needed for the plots, see below*) ThetaTable = {};(*values needed for recursive calculations, saves time*) For[t = 0, t =E2=89=A4 1000, t++,(*begin loop for the dynamical system consiting of t = 1000 Periods*) Clear[ThetaMatrix]; ThetaMatrix = Table[theta[i, k], {i, 1, 2}, {k, 1, 2}]; AppendTo[ThetaTable, ThetaMatrix]; AppendTo[RelWel, Rho[[1, 1]]*Wel[[1]]/Total[Rho.Wel]]; Clear[Wel]; x[t] = Random[];(*random pertubations are drawn from the above matrix*) If[x[t] =E2=89=A4 1/3, divi = Div[[1, All]], If[1/3 < x[t] =E2=89= =A4 2/3, divi = Div[[2, All]], divi = Div[[3, All]]]]; Wel = N[Inverse[IdentityMatrix[2] - ThetaMatrix.LamRho, Method -> DivisionFreeRowReduction].ThetaMatrix.divi, 1000];(*the inverse needed for t + 1*) ];(*end loop dynamical system*) Do[AppendTo[Value[t], RelWel[[t]]], {t, 1, 1000}]; ListPlot[RelWel, PlotStyle -> {Hue[.8]}, PlotRange -> {0, 1}, ImageSize -> 350] Clear[RelWel, ThetaTable]; ] (*end loop simulations*) ]