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

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.

• The above code is not fine since it evaluated with message NDSolve::ivcon: The given initial conditions were not consistent with the differential-algebraic equations. NDSolve will attempt to correct the values. Feb 6, 2022 at 12:33
• @AlexTrounev Thanks. I think my confusion has nothing to do with NDSolve, eventhoug Mathematica gives such a nice warning. solGolbal  and solModule use the same solution output by the NDSolve. In my option, even if this output is not the accurate solution to my ODEs, solGolbal  and solModule should still give the same reuslt. Feb 6, 2022 at 12:51
• The sample is not properly simplified, but the underlying problem is interesting. It can be boiled down to the following: Module[{rule = a -> b}, expr = Sin[t]; f[t_] = D[expr, t] /. rule] I myself didn't know that renaming will happen in this case. Since Module and renaming have been discussed quite a bit in this site, I won't be surprised if it's still a duplicate, but I don't have time to search for it at the moment. Feb 6, 2022 at 13:39
• @xzczd Yes, yes! This is the exact what I am looking for. I will search Module realted posts. Thanks! Feb 6, 2022 at 13:45
• A further simplification: Module[{rule = a -> b}, f[t_] = t /. rule] gives t\$. Feb 6, 2022 at 13:49

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.

• Thanks! I think I need more time to digest your helpful answer. Feb 6, 2022 at 14:11

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]


• Thanks. But this is not what I am asking for. I am trying to find the source of difference as shown in my figure. Feb 6, 2022 at 13:11
• dvarsol = D[varsol, t]; ` is bad programing style. Proposed a better way. Feb 6, 2022 at 14:06
• Yes, indeed, it is really bad coding style. Thanks! Feb 6, 2022 at 14:13