# Discrepancies in NDSolve Solution Between 4D and its Reduced 3D Counterpart

I have two systems, a 4D system and a 3D system. The 3D system is a reduction of the 4D system in such a way that:

(*
y1[t] = x2[t] - x1[t]
Thus:
-y1[t] = x1[t] - x2[t]
y1'[t] = x2'[t] - x1'[t]

y2[t] = x3[t] - x1[t]
Thus:
-y2[t] = x1[t] - x3[t]
y2'[t] = x3'[t] - x1'[t]

y3[t] = x4[t] - x1[t]
Thus:
-y1[t] = x1[t] - x4[t]
y1'[t] = x4'[t] - x1'[t]
*)

w = 1; k = 1/5; r = 1/4; a = 5/4; n = 4;

eqs = {
x1'[t]==w+(k/n)*(Sin[x2[t]-x1[t]+a]-r*Sin[2*(x2[t]-x1[t])] + Sin[x3[t]-x1[t]+a]-r*Sin[2*(x3[t]-x1[t])] + Sin[x4[t]-x1[t]+a]-r*Sin[2*(x4[t]-x1[t])]),

x2'[t]==w+(k/n)*(Sin[x1[t]-x2[t]+a]-r*Sin[2*(x1[t]-x2[t])] + Sin[x3[t]-x2[t]+a]-r*Sin[2*(x3[t]-x2[t])] + Sin[x4[t]-x2[t]+a]-r*Sin[2*(x4[t]-x2[t])]),

x3'[t]==w+(k/n)*(Sin[x1[t]-x3[t]+a]-r*Sin[2*(x1[t]-x3[t])] + Sin[x2[t]-x3[t]+a]-r*Sin[2*(x2[t]-x3[t])] + Sin[x4[t]-x3[t]+a]-r*Sin[2*(x4[t]-x3[t])]),

x4'[t]==w+(k/n)*(Sin[x1[t]-x4[t]+a]-r*Sin[2*(x1[t]-x4[t])] + Sin[x2[t]-x4[t]+a]-r*Sin[2*(x2[t]-x4[t])] + Sin[x3[t]-x4[t]+a]-r*Sin[2*(x3[t]-x4[t])])
};

eqsPD = {
y1'[t]==(k/n)*((Sin[-y1[t]+a]-r*Sin[-2*y1[t]]+Sin[y2[t]-y1[t]+a]-r*Sin[2*(y2[t]-y1[t])]+Sin[y3[t]-y1[t]+a]-r*Sin[2*(y3[t]-y1[t])]) - (Sin[y1[t]+a]-r*Sin[2*y1[t]]+Sin[y2[t]+a]-r*Sin[2*y2[t]]+Sin[y3[t]+a]-r*Sin[2*y3[t]])),

y2'[t]==(k/n)*((Sin[-y2[t]+a]-r*Sin[-2*y2[t]]+Sin[y1[t]-y2[t]+a]-r*Sin[2*(y1[t]-y2[t])]+Sin[y3[t]-y2[t]+a]-r*Sin[2*(y3[t]-y2[t])]) - (Sin[y1[t]+a]-r*Sin[2*y1[t]]+Sin[y2[t]+a]-r*Sin[2*y2[t]]+Sin[y3[t]+a]-r*Sin[2*y3[t]])),

y3'[t]==(k/n)*((Sin[-y3[t]+a]-r*Sin[-2*y3[t]]+Sin[y1[t]-y3[t]+a]-r*Sin[2*(y1[t]-y3[t])]+Sin[y2[t]-y3[t]+a]-r*Sin[2*(y2[t]-y3[t])]) - (Sin[y1[t]+a]-r*Sin[2*y1[t]]+Sin[y2[t]+a]-r*Sin[2*y2[t]]+Sin[y3[t]+a]-r*Sin[2*y3[t]]))
};


That is, I expect the solution for y1 to look like that of x2-x1, y2 to look like that of x3-x1, and y3 to look like that of x4-x1. However, solving them as below and plotting the solution gives:

ics={x1[0]==0, x2[0]==0.1, x3[0]==0.2, x4[0]==0.3};
icsPD={y1[0]==0.1, y2[0]==0.2, y3[0]==0.3};
eqsics=Join[eqs, ics];
eqsicsPD=Join[eqsPD, icsPD];
SolutionValue[t_]=NDSolveValue[eqsics,{x1[t],x2[t],x3[t],x4[t]},{t,0,15000}];
SolutionValuePD[t_]=NDSolveValue[eqsicsPD,{y1[t],y2[t],y3[t]},{t,0,15000}];


Where I have plotted x4-x1 in blue and y3 in dashed orange. I have tried increasing the AccuracyGoal to infinity so that I can work with the WorkingPrecision and PrecisionGoal instead, but doing this seems to always return that A component of the solution at t==0 is essentially zero and the tolerance specified by the AccuracyGoal option could not be achieved. Attempting to set WorkingPrecision returns that The precision of the differential equation is less than the WorkingPrecision for any value I seem to try.

In general, inconsistencies between the two systems seem to occur for values of a between 1.1 and 1.5, but otherwise they seem to be identical. I know that there are some SNIC or Hopf bifurcations in this parameter range that I imagine are related to why the system is so difficult to integrate for these values. Are there other approaches to solving this system more accurately that I am missing? Changing MaxStepSize to something smaller like 0.0001 has helped, however the computation starts to become rather expensive. Is continually decreasing the step size the best option or are there better approaches? And why / how is this occurring? Thank you for your time.