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.

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. — Alex Trounev, Feb 06 '22 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. — xinxin guo, Feb 06 '22 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. — xzczd, Feb 06 '22 at 13:39
@xzczd Yes, yes! This is the exact what I am looking for. I will search Module realted posts. Thanks! — xinxin guo, Feb 06 '22 at 13:45
A further simplification: Module[{rule = a -> b}, f[t_] = t /. rule] gives t$. — xzczd, Feb 06 '22 at 13:49
@xzczd Thanks for your constructive comments as always! The lesson I learned from this post is that : be becareful to use things like D[expr,t]. — xinxin guo, Feb 06 '22 at 14:15
I think that this is a problem of local and global variables. dpdxModule[t] not consists of -1000 D[U, t] while dpdxGlobal[t] consists it. We can improve this code by removing [t_] and [t] from {dpdxModule], dpdxGlobal}. — Alex Trounev, Feb 06 '22 at 14:17
I believe the following is a minimal example: Module[{x}, f[t_] = t x] This behavior is already noticed at least in https://mathematica.stackexchange.com/q/119915/1871 — xzczd, Feb 06 '22 at 15:16

score 3 · Answer 1 · edited Feb 06 '22 at 14:20

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. — xinxin guo, Feb 06 '22 at 13:11
dvarsol = D[varsol, t]; is bad programing style. Proposed a better way. — Akku14, Feb 06 '22 at 14:06

score 3 · Accepted Answer · answered Feb 06 '22 at 13:56

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. — xinxin guo, Feb 06 '22 at 14:11

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

2 Answers2