1

I am simulating the dynamics of a spring mass system impact an obstacle. I will not go into the details of the numerical scheme, but here are the general parameters and functions:

n = 10;
m = 1/n*IdentityMatrix[n];
k = n*SparseArray[{{i_, j_} /; Abs[i - j] == 1 -> -1, {i_, i_} /; 
        i != n -> 2, {n, n} -> 1}, {n, n}] // Normal;
c = 0.02*k;
F0 = 10;
e = 1;
d = 1;
T = 3.4;
h = 0.002;
theta = .5;
mhat = m + h*theta*c + h^2*theta^2*k;
mhatinv = Inverse[mhat];
Fext[i_] := {F0*Cos[2 Pi/T*h*i]}~Join~ConstantArray[0, n - 1];
prox[x_, b_] := Piecewise[{{-b, x < d}, {Min[0, -b], x >= d}}]

Then I define the following 3 sequences:

Clear[ulist, vlist, vfree];

ulist[0] = ConstantArray[0, n];
vlist[0] = ConstantArray[0, n];
vfree[0] = vlist[0];

vfree[i_] := vlist[i - 1] + mhatinv.(-h*c.vlist[i - 1] - h*k.ulist[i - 1] - 
     h^2*theta*k.vlist[i - 1] + h*(theta*Fext[i] + (1 - theta)*Fext[i - 1]))
ulist[i_] := ulist[i - 1] + h*(theta*vlist[i] + (1 - theta)*vlist[i - 1])
vlist[i_] := vfree[i][[;; -2]]~Join~{-e*vlist[i - 1][[-1]] + (1 + e)*
     prox[ulist[i - 1][[n]], -vfree[i][[n]]]}

Then it is possible to simulate the time-evolution of the mechanical system, and plot e.g. the trajectory of the contacting mass:

Table[ulist[i] = ulist[i]; vfree[i] = vfree[i]; vlist[i] = vlist[i], {i, 3000}];
Show[
   ListPlot[Table[{i*h, ulist[i][[n]]}, {i, 3000}]], 
   Plot[1, {x, 0, 10}, PlotStyle -> Red], 
   AxesLabel -> {"Time", "Position"}]

enter image description here

It works pretty well, however I am pretty sure this code can be improved a lot in terms of computational efficiency. I am thinking of using Compile, but it is not clear to me what to include (e.g., should the functions prox and Fext be defined in Compile? Should Compile be called for each iteration, or should the iteration process be included in it?).

Any hints would be appreciated.

anderstood
  • 14,301
  • 2
  • 29
  • 80

1 Answers1

4

I would start by memoizing these functions:

mem : vfree[i_] := mem =
  vlist[i - 1] + 
   mhatinv.(-h*c.vlist[i - 1] - h*k.ulist[i - 1] - h^2*theta*k.vlist[i - 1] + 
      h*(theta*Fext[i] + (1 - theta)*Fext[i - 1]))

mem : ulist[i_] := mem =
  ulist[i - 1] + h*(theta*vlist[i] + (1 - theta)*vlist[i - 1])

mem : vlist[i_] := mem =
  vfree[i][[;; -2]]~
   Join~{-e*vlist[i - 1][[-1]] + (1 + e)*prox[ulist[i - 1][[n]], -vfree[i][[n]]]}

Then eliminate this code as it becomes redundant:

Table[ulist[i] = ulist[i]; vfree[i] = vfree[i]; vlist[i] = vlist[i], {i, 30000}];

For me this takes total plotting time from 4.5 seconds down to 0.4 second.

Mr.Wizard
  • 271,378
  • 34
  • 587
  • 1,371
  • What is the function of mem here? Is it just a shorthand for repeating the function name in f[x_] := f[x] = rhs? – aardvark2012 Aug 15 '17 at 02:28
  • 2
    @aardvark2012 Yes, just personal convention; see the bottom of (2676) for my reasoning. – Mr.Wizard Aug 15 '17 at 02:29
  • Nice. I read the nice post you mention in the comment above, but don't understand why it is so much faster than when saving f[x]=f[x] afterwards. Could you please elaborate a bit on that? Also, do you think it could be still possible to increase the speed significantly? – anderstood Aug 15 '17 at 15:07
  • 1
    @anderstood (1) I think because only values used are computed. I didn't check to see if that was equivalent in this case, and looking now I see that the plot is only to 3,000 rather than the 30,000 used in the Table. That by itself would account for the timings I observed. In some cases though the order of computation is important. (2) I hope so. As I said I would start by ... but I didn't have time to look for other optimizations. I did note that the output of ulist and vlist are already packed arrays which is one sign of optimized code. – Mr.Wizard Aug 15 '17 at 17:47