# Algorithm to approximate non-linear equation system solution

I'm looking for an algorithm to approximate the solution of the following equation system:

The equations have to be solved on an embedded system, in C++.

Background:

• We measure the 2 variables X_m and Y_m, so they are known
• We want to compute the real values: X_r and Y_r
• X and Y are real numbers
• We measure the functions f_xy and f_yx during calibration. We have maximal 18 points of each function.
• It's possible to store the functions as a look-up table
• I tried to approximate the functions with 2nd order polynomials and compute the solution, but it was not accurate enough, because of the fitting error.

I am looking for an algorithm to approximate the results in an embedded system in C++, but I don't even know what to search for. I found some papers on the theory link, but I think there must be an easier way to do it in my case.

Also: how can I determine during calibration, whether the functions can be solved with the algorithm?

The normal way is by Newton's iterations, starting from the initial approximation `(Xm, Ym)` [assuming that the `f` are mere corrections]. Due to the particular shape of the equations, you can reduce to twice a single equation in a single unknown.

``````Xr = Xm - Fyx(Ym - Fxy(Xr))
Yr = Ym - Fxy(Xm - Fyx(Yr))
``````

``````Xr <-- Xr - (Xm - Fyx(Ym - Fxy(Xr))) / (1 + Fxy'(Ym - Fxy(Xr)).Fxy'(Xr))
Yr <-- Yr - (Ym - Fxy(Xm - Fyx(Yr))) / (1 + Fyx'(Xm - Fyx(Yr)).Fyx'(Yr))
``````

So you should tabulate the derivatives of the `f` as well, though accuracy is not so critical than for the computation of the `f` themselves.

If the calibration points aren't too noisy, I would recommend cubic spline interpolation, for which you can precompute all coefficients. At the same time these coefficients allow you to estimate the derivative (as the corresponding quadratic interpolant, which is continuous).

In principle (unless the points are uniformly spaced), you need to perform a dichotomic search to determine the interval in which the argument lies. But here you will evaluate the functions at nearby values, so that a linear search from the previous location should be better.

A different way to address the problem is by considering the bivariate solution surfaces `Xr = G(Xm, Ym)` and `Yr = G(Xm, Ym)` that you compute on a grid of points. If the surfaces are smooth enough, you can use a coarse grid.

So by any method (such as the one above), you precompute the solutions at each grid node, as well as the coefficients of some interpolant in the `X` and `Y` directions. I recommend a cubic spline, again.

Now to interpolate inside a grid cell, you combine the two univarite interpolants to a bivariate one by means of the Coons formula https://en.wikipedia.org/wiki/Coons_patch.

Fitting a second-order polynomial through `f_xy`? That's generally not viable. The go-to solution would be Runga-Kutta interpolation. You pick two known values left and two to the right of your argument, with weights 1,2,2,1. This gets you an estimate `d(f_xy)/dx` which you can then use for interpolation.