In the context of Quantum Computing, I am investigating curves that illustrate an error propagation during qubit rotations on a Blochsphere (I asked a related question in the same context here and received an useful answer). The plot below depicts two of such curves:

enter image description here

Both curves show an angular deviation in azimuth and elevation (that is the reason that we have two curves). The code is not that complicated, however it utilizes function nesting:

qubitmatrixToCartesian[Mq_] := (
   q1 = Re[(Mq[[1, 2]] + Mq[[2, 1]])/2];
   q2 = Re[(Mq[[2, 1]] - Mq[[1, 2]])/(2*I)];
   q3 = Re[Mq[[1, 1]]];
   Return[{q1, q2, q3}];

rnSU2euler[vec_, rx_, ry_, rz_] := (
   sphericalVec = ToSphericalCoordinates[vec];
   \[Theta] = sphericalVec[[2]];
   \[Phi] = sphericalVec[[3]];
   sx = PauliMatrix[1];
   sy = PauliMatrix[2];
   sz = PauliMatrix[3];
   Mq = Sin[\[Theta]]*Cos[\[Phi]]*sx + Sin[\[Theta]]*Sin[\[Phi]]*sy + 
   Un = {{Exp[-I*(rx + rz)/2]*Cos[ry/2], -Exp[-I*(rx - rz)/2]*
       Sin[ry/2]}, {Exp[I*(rx - rz)/2]*Sin[ry/2], 
      Exp[I*(rx + rz)/2]*Cos[ry/2]}};
   Return [Un . Mq . ConjugateTranspose[Un]];

rotateVector[vec_, rx_, ry_, rz_] := (
   MqRotated = rnSU2euler[vec, rx, ry, rz];

errorPropagation[n_, vec_, vecError_, rx_, ry_, rz_] := (
   v = NestList[rotateVector[#, rx, ry, rz] &, N@vec, n];
   vErr = NestList[rotateVector[#, rx, ry, rz] &, N@vecError, n];
   polar = Map[ToSphericalCoordinates, v];
   polarErr = Map[ToSphericalCoordinates, vErr];
   {polarErr[[All, 3]] - polar[[All, 3]], 
    polarErr[[All, 2]] - polar[[All, 2]]}

 errorPropagation[200, {1, 0, 0}, 
  N[rotateVector[{1, 0, 0}, 0, 0.2, 0]], Pi/100, Pi/100, Pi/100]]

Is it possible for Mathematica to determine the curve equation of these two curves? I know that Mathematica is quite powerful, as even it can solve recursive functions (via RSolve). I would appreciate any advice on how to examine the curves. Even if it is not possible to determine the curve equation, I would be grateful for any hints, for example, to allow determining other characteristics such as periodicity, etc.

  • 1
    $\begingroup$ The first curve is discontinuous. What sort of equation do you think about. Piecewise? $\endgroup$ Jun 25, 2022 at 17:55
  • $\begingroup$ @DanielHuber: Thank you for your hint. Yes - for me even a Piecewise function would be fine as well. Note: Possibly both curves are continuous - here is a full plot. $\endgroup$ Jun 25, 2022 at 18:13
  • 1
    $\begingroup$ The simplest you could do is linear interpolation. The draw back is you will have kinks. If you do not need the explicit formula, you may use "Interpolation" with piecewise polynomial or spline interpolation.. $\endgroup$ Jun 25, 2022 at 19:27

1 Answer 1


You would have to come up with a model in order to fit those curves. No automated system can do that for you. However, that may not be the most useful approach. For instance some properties of the data (e.g. periodicity) can be extracted from the raw data by Fourier transform using e.g. Periodogram.

Let's generate some data using your function definitions:

data = errorPropagation[
         30000, {1, 0, 0}, 
         N[rotateVector[{1, 0, 0}, 0, 0.2, 0]],
         Pi/100, Pi/100, Pi/100

  data[[All, ;; 1000]], Joined -> True, 
   PlotLegends -> {"first", "second"}

two data sets in blue (first) and yellow (second)

The second curve in this set, in yellow in the plot above, seems easier to tackle:

  data[[2, ;; 500]],
  PlotRange -> All, PlotStyle -> ColorData[97][2]

second curve only, zoomed in

Its periodogram is:

  PlotRange -> {{1*^-3, 0.11}, All}, 
  Frame -> True, Axes -> False,
  PlotStyle -> ColorData[97][2]

periodogram of yellow curve

The evenly spaces peaks hint at harmonics of the primary frequency, at around 0.01 Hz. We are mostly interested in the period, i.e. the inverse of the frequency, which in your case we can understand as "iterations" or "cycles", similar to the first argument to your errorPropagation function.

To estimate the position of that first peak, you could use FindPeaks, but in this case I will just zooming in to that region of the plot and read it off the chart using "Get Coordinates...":

  data[[2]], PlotRange -> {{0.0108, 0.0115}, {-1, 8}},
  ScalingFunctions -> None,
  Frame -> True, Axes -> False, PlotStyle -> ColorData[97][2]

zoomed-in periodogram

My best estimate of the peak's location is {{0.011166709533541064, 7.1772552649066235}}. We don't care about the intensity; only about the period, which is 1/0.011166709533541064 == 89.5519. So the period seems to be 89.5 points/iterations or thereabouts.

We can visually check that on a plot of the raw data by inserting some GridLines with that spacing:

  data[[2, ;; 1000]],
  GridLines -> {Range[0, 1000, 89.55], None},
  GridLinesStyle -> Black, PlotStyle -> ColorData[97][2]

raw data with gridlines indicating periodicity

Hopefully a similar approach can be used for the data[[1]] dataset as well.

  • 1
    $\begingroup$ Thank you for pointing me to the Periodogram. This is interesting and I wasn't aware of this feature. $\endgroup$ Jun 25, 2022 at 19:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.