Replacement with equivalent local variable and global variables give apparently different results using the same NDSolve output

Question

I am solving a system ODEs which Mathematica quickly gives me a result varsol and dvarsol. Up to this point, everything is good.

However, when I substitute these two solution into my final expression (-1000 D[U, t] + S)[[5]], I found a strange thing.

I tried two ways of substituting which should have be equivalent.

(1) The fist way is that I collected varsol and dvarsol in a global variable solGolbal then do the substituting;

(2)The second way is that I collected varsol and dvarsol in a local variable solModule and then do the substituting.

The results of these two ways I think should be the same, however they showd apparent difference (i.e., the global way has one more term than the local way, as shown in the red box in the figure below).

Remove["Global`*"] // Quiet;
sys = {s[1][t] + s[1][t]^2/40000 + 1/100 Derivative[1][s[1]][t] == 
    1/50 \[Tau][1][
      t] (256000000/9 \[CurlyPhi][3][t] - 
       5120/9 (-9 \[Pi] Sin[64 \[Pi] t] + 15625 \[CurlyPhi][3][t]) - 
       64000000/9 \[CurlyPhi][4][t]), 
   s[2][t] + s[2][t]^2/40000 + 1/100 Derivative[1][s[2]][t] == 
    1/50 \[Tau][2][
      t] (64000000/9 \[CurlyPhi][3][t] - 
       2048/9 (-9 \[Pi] Sin[64 \[Pi] t] + 15625 \[CurlyPhi][3][t])), 
   s[3][t] + s[3][t]^2/40000 + 1/100 Derivative[1][s[3]][t] == 
    1/50 \[Tau][3][
      t] (-(128000000/9) \[CurlyPhi][3][t] + 
       1024/9 (-9 \[Pi] Sin[64 \[Pi] t] + 15625 \[CurlyPhi][3][t]) + 
       64000000/9 \[CurlyPhi][4][t]), 
   s[4][t] + s[4][t]^2/40000 + 1/100 Derivative[1][s[4]][t] == 
    1/50 \[Tau][4][
      t] (64000000/9 \[CurlyPhi][3][t] - 
       128000000/9 \[CurlyPhi][4][t] + 64000000/9 \[CurlyPhi][5][t]), 
   s[5][t] + s[5][t]^2/40000 + 1/100 Derivative[1][s[5]][t] == 
    1/50 \[Tau][5][
      t] (64000000/9 \[CurlyPhi][4][t] - 
       128000000/9 \[CurlyPhi][5][t] + 64000000/9 \[CurlyPhi][6][t]), 
   s[6][t] + s[6][t]^2/40000 + 1/100 Derivative[1][s[6]][t] == 
    1/50 \[Tau][6][
      t] (64000000/9 \[CurlyPhi][5][t] - 
       128000000/9 \[CurlyPhi][6][t] + 64000000/9 \[CurlyPhi][7][t]), 
   s[7][t] + s[7][t]^2/40000 + 1/100 Derivative[1][s[7]][t] == 
    1/50 \[Tau][7][
      t] (64000000/9 \[CurlyPhi][6][t] + 
       64000000/
        9 ((4131 \[Pi] Sin[64 \[Pi] t])/734375 + 
          1/4 \[CurlyPhi][7][t]) - 128000000/9 \[CurlyPhi][7][t]), 
   s[8][t] + s[8][t]^2/40000 + 1/100 Derivative[1][s[8]][t] == 
    1/50 \[Tau][8][
      t] (2506752/47 \[Pi] Sin[64 \[Pi] t] - 
       128000000/
        9 ((4131 \[Pi] Sin[64 \[Pi] t])/734375 + 
          1/4 \[CurlyPhi][7][t]) + 64000000/9 \[CurlyPhi][7][t]), 
   s[9][t] + s[9][t]^2/40000 + 1/100 Derivative[1][s[9]][t] == 
    1/50 \[Tau][9][
      t] (5013504/47 \[Pi] Sin[64 \[Pi] t] - 
       64000000/9 \[CurlyPhi][6][t] - 
       320000000/
        9 ((4131 \[Pi] Sin[64 \[Pi] t])/734375 + 
          1/4 \[CurlyPhi][7][t]) + 
       256000000/9 \[CurlyPhi][7][t]), \[Tau][1][t] + (
     s[1][t] \[Tau][1][t])/40000 + 
     1/100 Derivative[1][\[Tau][1]][t] == 
    100 (256000000/9 \[CurlyPhi][3][t] - 
       5120/9 (-9 \[Pi] Sin[64 \[Pi] t] + 15625 \[CurlyPhi][3][t]) - 
       64000000/9 \[CurlyPhi][4][t]), \[Tau][2][t] + (
     s[2][t] \[Tau][2][t])/40000 + 
     1/100 Derivative[1][\[Tau][2]][t] == 
    100 (64000000/9 \[CurlyPhi][3][t] - 
       2048/
        9 (-9 \[Pi] Sin[64 \[Pi] t] + 
          15625 \[CurlyPhi][3][t])), \[Tau][3][t] + (
     s[3][t] \[Tau][3][t])/40000 + 
     1/100 Derivative[1][\[Tau][3]][t] == 
    100 (-(128000000/9) \[CurlyPhi][3][t] + 
       1024/9 (-9 \[Pi] Sin[64 \[Pi] t] + 15625 \[CurlyPhi][3][t]) + 
       64000000/9 \[CurlyPhi][4][t]), \[Tau][4][t] + (
     s[4][t] \[Tau][4][t])/40000 + 
     1/100 Derivative[1][\[Tau][4]][t] == 
    100 (64000000/9 \[CurlyPhi][3][t] - 
       128000000/9 \[CurlyPhi][4][t] + 
       64000000/9 \[CurlyPhi][5][t]), \[Tau][5][t] + (
     s[5][t] \[Tau][5][t])/40000 + 
     1/100 Derivative[1][\[Tau][5]][t] == 
    100 (64000000/9 \[CurlyPhi][4][t] - 
       128000000/9 \[CurlyPhi][5][t] + 
       64000000/9 \[CurlyPhi][6][t]), \[Tau][6][t] + (
     s[6][t] \[Tau][6][t])/40000 + 
     1/100 Derivative[1][\[Tau][6]][t] == 
    100 (64000000/9 \[CurlyPhi][5][t] - 
       128000000/9 \[CurlyPhi][6][t] + 
       64000000/9 \[CurlyPhi][7][t]), \[Tau][7][t] + (
     s[7][t] \[Tau][7][t])/40000 + 
     1/100 Derivative[1][\[Tau][7]][t] == 
    100 (64000000/9 \[CurlyPhi][6][t] + 
       64000000/
        9 ((4131 \[Pi] Sin[64 \[Pi] t])/734375 + 
          1/4 \[CurlyPhi][7][t]) - 
       128000000/9 \[CurlyPhi][7][t]), \[Tau][8][t] + (
     s[8][t] \[Tau][8][t])/40000 + 
     1/100 Derivative[1][\[Tau][8]][t] == 
    100 (2506752/47 \[Pi] Sin[64 \[Pi] t] - 
       128000000/
        9 ((4131 \[Pi] Sin[64 \[Pi] t])/734375 + 
          1/4 \[CurlyPhi][7][t]) + 
       64000000/9 \[CurlyPhi][7][t]), \[Tau][9][t] + (
     s[9][t] \[Tau][9][t])/40000 + 
     1/100 Derivative[1][\[Tau][9]][t] == 
    100 (5013504/47 \[Pi] Sin[64 \[Pi] t] - 
       64000000/9 \[CurlyPhi][6][t] - 
       320000000/
        9 ((4131 \[Pi] Sin[64 \[Pi] t])/734375 + 
          1/4 \[CurlyPhi][7][t]) + 256000000/9 \[CurlyPhi][7][t]), 
   128000000/9 Derivative[1][\[CurlyPhi][3]][t] - 
     1024/
      9 (-576 \[Pi]^2 Cos[64 \[Pi] t] + 
        15625 Derivative[1][\[CurlyPhi][3]][t]) - 
     64000000/
      9 Derivative[1][\[CurlyPhi][4]][t] == (-(64000000/9) \[Tau][2][
       t] + 128000000/9 \[Tau][3][t] - 64000000/9 \[Tau][4][t])/
    1000, -(64000000/9) Derivative[1][\[CurlyPhi][3]][t] + 
     128000000/9 Derivative[1][\[CurlyPhi][4]][t] - 
     64000000/
      9 Derivative[1][\[CurlyPhi][5]][t] == (-(64000000/9) \[Tau][3][
       t] + 128000000/9 \[Tau][4][t] - 64000000/9 \[Tau][5][t])/
    1000, -(64000000/9) Derivative[1][\[CurlyPhi][4]][t] + 
     128000000/9 Derivative[1][\[CurlyPhi][5]][t] - 
     64000000/
      9 Derivative[1][\[CurlyPhi][6]][t] == (-(64000000/9) \[Tau][4][
       t] + 128000000/9 \[Tau][5][t] - 64000000/9 \[Tau][6][t])/
    1000, -(64000000/9) Derivative[1][\[CurlyPhi][5]][t] + 
     128000000/9 Derivative[1][\[CurlyPhi][6]][t] - 
     64000000/
      9 Derivative[1][\[CurlyPhi][7]][t] == (-(64000000/9) \[Tau][5][
       t] + 128000000/9 \[Tau][6][t] - 64000000/9 \[Tau][7][t])/
    1000, -(64000000/9) Derivative[1][\[CurlyPhi][6]][t] - 
     64000000/
      9 ((264384 \[Pi]^2 Cos[64 \[Pi] t])/734375 + 
        1/4 Derivative[1][\[CurlyPhi][7]][t]) + 
     128000000/
      9 Derivative[1][\[CurlyPhi][7]][t] == (-(64000000/9) \[Tau][6][
       t] + 128000000/9 \[Tau][7][t] - 64000000/9 \[Tau][8][t])/1000, 
   s[1][0] == 0, s[2][0] == 0, s[3][0] == 0, s[4][0] == 0, 
   s[5][0] == 0, s[6][0] == 0, s[7][0] == 0, s[8][0] == 0, 
   s[9][0] == 0, \[Tau][1][0] == 0, \[Tau][2][0] == 0, \[Tau][3][0] ==
     0, \[Tau][4][0] == 0, \[Tau][5][0] == 0, \[Tau][6][0] == 
    0, \[Tau][7][0] == 0, \[Tau][8][0] == 0, \[Tau][9][0] == 
    0, \[CurlyPhi][3][0] == 0, \[CurlyPhi][4][0] == 
    0, \[CurlyPhi][5][0] == 0, \[CurlyPhi][6][0] == 
    0, \[CurlyPhi][7][0] == 0, Derivative[1][s[1]][0] == 1, 
   Derivative[1][s[2]][0] == 1, Derivative[1][s[3]][0] == 1, 
   Derivative[1][s[4]][0] == 1, Derivative[1][s[5]][0] == 1, 
   Derivative[1][s[6]][0] == 1, Derivative[1][s[7]][0] == 1, 
   Derivative[1][s[8]][0] == 1, Derivative[1][s[9]][0] == 1, 
   Derivative[1][\[Tau][1]][0] == 1/100000000, 
   Derivative[1][\[Tau][2]][0] == 1/100000000, 
   Derivative[1][\[Tau][3]][0] == 1/100000000, 
   Derivative[1][\[Tau][4]][0] == 1/100000000, 
   Derivative[1][\[Tau][5]][0] == 1/100000000, 
   Derivative[1][\[Tau][6]][0] == 1/100000000, 
   Derivative[1][\[Tau][7]][0] == 1/100000000, 
   Derivative[1][\[Tau][8]][0] == 1/100000000, 
   Derivative[1][\[Tau][9]][0] == 1/100000000, 
   Derivative[1][\[CurlyPhi][3]][0] == 1, 
   Derivative[1][\[CurlyPhi][4]][0] == 1, 
   Derivative[1][\[CurlyPhi][5]][0] == 1, 
   Derivative[1][\[CurlyPhi][6]][0] == 1, 
   Derivative[1][\[CurlyPhi][7]][0] == 1};

var = {\[CurlyPhi][3][t], \[CurlyPhi][4][t], \[CurlyPhi][5][
    t], \[CurlyPhi][6][t], \[CurlyPhi][7][t], s[1][t], s[2][t], 
   s[3][t], s[4][t], s[5][t], s[6][t], s[7][t], s[8][t], 
   s[9][t], \[Tau][1][t], \[Tau][2][t], \[Tau][3][t], \[Tau][4][
    t], \[Tau][5][t], \[Tau][6][t], \[Tau][7][t], \[Tau][8][
    t], \[Tau][9][t]};

{U, S} = {{-(96/125) \[Pi] Sin[64 \[Pi] t], 4000/3 \[CurlyPhi][3][t], 
    8/375 (9 \[Pi] Sin[64 \[Pi] t] - 15625 \[CurlyPhi][3][t] + 
       62500 \[CurlyPhi][4][t]), -(4000/
      3) (\[CurlyPhi][3][t] - \[CurlyPhi][5][t]), -(4000/
      3) (\[CurlyPhi][4][t] - \[CurlyPhi][6][t]), -(4000/
      3) (\[CurlyPhi][5][t] - \[CurlyPhi][7][t]), (
     44064 \[Pi] Sin[64 \[Pi] t])/5875 - 4000/3 \[CurlyPhi][6][t] + 
     1000/3 \[CurlyPhi][7][t], (58752 \[Pi] Sin[64 \[Pi] t])/5875 - 
     4000/3 \[CurlyPhi][7][t], 
    0}, {-4000 \[Tau][1][t] + 16000/3 \[Tau][2][t] - 
     4000/3 \[Tau][3][t], -(4000/3) \[Tau][1][t] + 
     4000/3 \[Tau][3][t], -(4000/3) \[Tau][2][t] + 
     4000/3 \[Tau][4][t], -(4000/3) \[Tau][3][t] + 
     4000/3 \[Tau][5][t], -(4000/3) \[Tau][4][t] + 
     4000/3 \[Tau][6][t], -(4000/3) \[Tau][5][t] + 
     4000/3 \[Tau][7][t], -(4000/3) \[Tau][6][t] + 
     4000/3 \[Tau][8][t], -(4000/3) \[Tau][7][t] + 
     4000/3 \[Tau][9][t], 
    4000/3 \[Tau][7][t] - 16000/3 \[Tau][8][t] + 4000 \[Tau][9][t]}};

T = 0.03125;
varsol = 
  NDSolve[sys, var, {t, 0, 6 T}, 
    Method -> {"EquationSimplification" -> "Residual"}] // Flatten;
dvarsol = D[varsol, t];
vp[t_] = -(96/125) \[Pi] Sin[64 \[Pi] t];

------- The above code is fine. The problem occured in following lines. ------

---

compute[_] := Module[{solModule},
   (*local way*)
   solModule = {varsol, dvarsol} // Flatten; 
   dpdxModule[t_] = Evaluate[(-1000 D[U, t] + S)[[5]] /. solModule];
   (*global way*)
   solGolbal = {varsol, dvarsol} // Flatten;
   dpdxGlobal[t_] = Evaluate[(-1000 D[U, t] + S)[[5]] /. solGolbal];
   Print[solModule - solGolbal];
   {dpdxModule[t], dpdxGlobal[t]}];

F[t_] = compute[1];


F[t](*global way has one more term than local way*)

ParametricPlot[{{vp[t], F[t][[1]]}, {vp[t], F[t][[2]]}}, {t, 5 T, 
  6 T}, AspectRatio -> 1/GoldenRatio, 
 PlotLegends -> {"use module solution", "use global solution"}]

The problem confusing me is :

Why did dpdxGlobal[t] have one more term than dpdxModule[t], as shown in the red box in my second figure?

Thanks.

Answer 1

This seems to be a problem from mixing local and global variable with the same name. The subsequent differences in numerical calculations seems to come from taking the differences of large numbers.

Add Print[{dpdxModule[0], dpdxGlobal[0]}]; to the code of compute. The output is:

As you can see, "t" in the first case is not replaced by zero. What is the reason. In my opinion the culprit is:

F[t_] = compute[1];

This is not very clean written. The "1" does nothing. And it is not clear if the "t" From F[t_] and the t in "compute" are the same. In "compute", which t is taken? the global, or the localized one from F[t_]? I think the global one is taken. Then we have, in U and S, a further "t", which is clearly global.

Further in "solModule" is "t" global or local? It think it is local. Therefore, it is not replaced in the Evaluate inside compute.

The resulting differences in numerical calculations seems to come from taking the differences of large numbers.

How to fix it? IMHO, you should clean up parameters in functions. Then it should work.

Answer 2

If you need derivatives of variables for later calculations, the best way is to solve for pure functions (Regard, how var is defined) , where you get derivative very simply as shown. Second, with first order equations, you may not give derivative boundary conditions.(I gidn't check your solution)

sys = {s[1][t] + s[1][t]^2/40000 + 1/100 Derivative[1][s[1]][t] == 
1/50 τ[1][
  t] (256000000/9 φ[3][t] - 
   5120/9 (-9 π Sin[64 π t] + 15625 φ[3][t]) - 
   64000000/9 φ[4][t]), 
s[2][t] + s[2][t]^2/40000 + 1/100 Derivative[1][s[2]][t] == 
1/50 τ[2][
  t] (64000000/9 φ[3][t] - 
   2048/9 (-9 π Sin[64 π t] + 15625 φ[3][t])), 
s[3][t] + s[3][t]^2/40000 + 1/100 Derivative[1][s[3]][t] == 
1/50 τ[3][
  t] (-(128000000/9) φ[3][t] + 
   1024/9 (-9 π Sin[64 π t] + 15625 φ[3][t]) + 
   64000000/9 φ[4][t]), 
s[4][t] + s[4][t]^2/40000 + 1/100 Derivative[1][s[4]][t] == 
1/50 τ[4][
  t] (64000000/9 φ[3][t] - 
   128000000/9 φ[4][t] + 64000000/9 φ[5][t]), 
s[5][t] + s[5][t]^2/40000 + 1/100 Derivative[1][s[5]][t] == 
1/50 τ[5][
  t] (64000000/9 φ[4][t] - 
   128000000/9 φ[5][t] + 64000000/9 φ[6][t]), 
s[6][t] + s[6][t]^2/40000 + 1/100 Derivative[1][s[6]][t] == 
1/50 τ[6][
  t] (64000000/9 φ[5][t] - 
   128000000/9 φ[6][t] + 64000000/9 φ[7][t]), 
s[7][t] + s[7][t]^2/40000 + 1/100 Derivative[1][s[7]][t] == 
1/50 τ[7][
  t] (64000000/9 φ[6][t] + 
   64000000/
     9 ((4131 π Sin[64 π t])/734375 + 
      1/4 φ[7][t]) - 128000000/9 φ[7][t]), 
s[8][t] + s[8][t]^2/40000 + 1/100 Derivative[1][s[8]][t] == 
1/50 τ[8][
  t] (2506752/47 π Sin[64 π t] - 
   128000000/
     9 ((4131 π Sin[64 π t])/734375 + 
      1/4 φ[7][t]) + 64000000/9 φ[7][t]), 
s[9][t] + s[9][t]^2/40000 + 1/100 Derivative[1][s[9]][t] == 
1/50 τ[9][
  t] (5013504/47 π Sin[64 π t] - 
   64000000/9 φ[6][t] - 
   320000000/
     9 ((4131 π Sin[64 π t])/734375 + 
      1/4 φ[7][t]) + 
   256000000/9 φ[7][t]), τ[1][
  t] + (s[1][t] τ[1][t])/40000 + 
 1/100 Derivative[1][τ[1]][t] == 
100 (256000000/9 φ[3][t] - 
   5120/9 (-9 π Sin[64 π t] + 15625 φ[3][t]) - 
   64000000/9 φ[4][t]), τ[2][
  t] + (s[2][t] τ[2][t])/40000 + 
 1/100 Derivative[1][τ[2]][t] == 
100 (64000000/9 φ[3][t] - 
   2048/9 (-9 π Sin[64 π t] + 
      15625 φ[3][t])), τ[3][
  t] + (s[3][t] τ[3][t])/40000 + 
 1/100 Derivative[1][τ[3]][t] == 
100 (-(128000000/9) φ[3][t] + 
   1024/9 (-9 π Sin[64 π t] + 15625 φ[3][t]) + 
   64000000/9 φ[4][t]), τ[4][
  t] + (s[4][t] τ[4][t])/40000 + 
 1/100 Derivative[1][τ[4]][t] == 
100 (64000000/9 φ[3][t] - 
   128000000/9 φ[4][t] + 
   64000000/9 φ[5][t]), τ[5][
  t] + (s[5][t] τ[5][t])/40000 + 
 1/100 Derivative[1][τ[5]][t] == 
100 (64000000/9 φ[4][t] - 
   128000000/9 φ[5][t] + 
   64000000/9 φ[6][t]), τ[6][
  t] + (s[6][t] τ[6][t])/40000 + 
 1/100 Derivative[1][τ[6]][t] == 
100 (64000000/9 φ[5][t] - 
   128000000/9 φ[6][t] + 
   64000000/9 φ[7][t]), τ[7][
  t] + (s[7][t] τ[7][t])/40000 + 
 1/100 Derivative[1][τ[7]][t] == 
100 (64000000/9 φ[6][t] + 
   64000000/
     9 ((4131 π Sin[64 π t])/734375 + 
      1/4 φ[7][t]) - 
   128000000/9 φ[7][t]), τ[8][
  t] + (s[8][t] τ[8][t])/40000 + 
 1/100 Derivative[1][τ[8]][t] == 
100 (2506752/47 π Sin[64 π t] - 
   128000000/
     9 ((4131 π Sin[64 π t])/734375 + 
      1/4 φ[7][t]) + 
   64000000/9 φ[7][t]), τ[9][
  t] + (s[9][t] τ[9][t])/40000 + 
 1/100 Derivative[1][τ[9]][t] == 
100 (5013504/47 π Sin[64 π t] - 
   64000000/9 φ[6][t] - 
   320000000/
     9 ((4131 π Sin[64 π t])/734375 + 
      1/4 φ[7][t]) + 256000000/9 φ[7][t]), 
128000000/9 Derivative[1][φ[3]][t] - 
 1024/9 (-576 π^2 Cos[64 π t] + 
    15625 Derivative[1][φ[3]][t]) - 
 64000000/9 Derivative[1][φ[4]][
   t] == (-(64000000/9) τ[2][t] + 128000000/9 τ[3][t] - 
   64000000/9 τ[4][t])/
 1000, -(64000000/9) Derivative[1][φ[3]][t] + 
 128000000/9 Derivative[1][φ[4]][t] - 
 64000000/9 Derivative[1][φ[5]][
   t] == (-(64000000/9) τ[3][t] + 128000000/9 τ[4][t] - 
   64000000/9 τ[5][t])/
 1000, -(64000000/9) Derivative[1][φ[4]][t] + 
 128000000/9 Derivative[1][φ[5]][t] - 
 64000000/9 Derivative[1][φ[6]][
   t] == (-(64000000/9) τ[4][t] + 128000000/9 τ[5][t] - 
   64000000/9 τ[6][t])/
 1000, -(64000000/9) Derivative[1][φ[5]][t] + 
 128000000/9 Derivative[1][φ[6]][t] - 
 64000000/9 Derivative[1][φ[7]][
   t] == (-(64000000/9) τ[5][t] + 128000000/9 τ[6][t] - 
   64000000/9 τ[7][t])/
 1000, -(64000000/9) Derivative[1][φ[6]][t] - 
 64000000/
   9 ((264384 π^2 Cos[64 π t])/734375 + 
    1/4 Derivative[1][φ[7]][t]) + 
 128000000/9 Derivative[1][φ[7]][
   t] == (-(64000000/9) τ[6][t] + 128000000/9 τ[7][t] - 
   64000000/9 τ[8][t])/1000, s[1][0] == 0, s[2][0] == 0, 
s[3][0] == 0, s[4][0] == 0, s[5][0] == 0, s[6][0] == 0, 
s[7][0] == 0, s[8][0] == 0, 
s[9][0] == 0, τ[1][0] == 0, τ[2][0] == 0, τ[3][0] ==
 0, τ[4][0] == 0, τ[5][0] == 0, τ[6][0] == 
 0, τ[7][0] == 0, τ[8][0] == 0, τ[9][0] == 
 0, φ[3][0] == 0, φ[4][0] == 
 0, φ[5][0] == 0, φ[6][0] == 
 0, φ[7][0] == 0};

.

var = {φ[3][t], φ[4][t], φ[5][
t], φ[6][t], φ[7][t], s[1][t], s[2][t], 
s[3][t], s[4][t], s[5][t], s[6][t], s[7][t], s[8][t], 
s[9][t], τ[1][t], τ[2][t], τ[3][t], τ[4][
t], τ[5][t], τ[6][t], τ[7][t], τ[8][
t], τ[9][t]} /. aa_[t] -> aa

({U, S} = {{-(96/125) π Sin[64 π t], 4000/3 φ[3][t],
  8/375 (9 π Sin[64 π t] - 15625 φ[3][t] + 
    62500 φ[4][t]), -(4000/3) (φ[3][
     t] - φ[5][t]), -(4000/3) (φ[4][
     t] - φ[6][t]), -(4000/3) (φ[5][
     t] - φ[7][t]), (44064 π Sin[64 π t])/
   5875 - 4000/3 φ[6][t] + 
  1000/3 φ[7][t], (58752 π Sin[64 π t])/5875 - 
  4000/3 φ[7][t], 
 0}, {-4000 τ[1][t] + 16000/3 τ[2][t] - 
  4000/3 τ[3][t], -(4000/3) τ[1][t] + 
  4000/3 τ[3][t], -(4000/3) τ[2][t] + 
  4000/3 τ[4][t], -(4000/3) τ[3][t] + 
  4000/3 τ[5][t], -(4000/3) τ[4][t] + 
  4000/3 τ[6][t], -(4000/3) τ[5][t] + 
  4000/3 τ[7][t], -(4000/3) τ[6][t] + 
  4000/3 τ[8][t], -(4000/3) τ[7][t] + 
  4000/3 τ[9][t], 
 4000/3 τ[7][t] - 16000/3 τ[8][t] + 
  4000 τ[9][t]}});

T = 0.03125 // Rationalize[#, 0] &;

varsol = NDSolve[sys, var, {t, 0, 6 T}] // Flatten;

.

vp[t_] = -(96/125) π Sin[64 π t];

ParametricPlot[
Evaluate[{vp[t], (-1000 D[U, t] + S)[[5]] /. varsol}], {t, 5 T, 6 T},
PlotStyle -> {Blue}, PlotRange -> Full, AspectRatio -> 1]