# How can I have Mathematica solve the curve equation when nested functions are involved?

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:

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:

ClearAll["Global*"];
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 +
Cos[\[Theta]]*sz;
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];
Return[qubitmatrixToCartesian[MqRotated]];
);

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]]}
);

ListLinePlot[
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.

• The first curve is discontinuous. What sort of equation do you think about. Piecewise? Jun 25, 2022 at 17:55
• @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. Jun 25, 2022 at 18:13
• 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.. Jun 25, 2022 at 19:27

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
];

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


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

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


Its periodogram is:

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


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...":

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


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:

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


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

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