I am coming from Stack overflow, where they closed my question because of being "mathemathical and nothing of programming", please, I know that in this forums usually mathemathical questions are more complex and elaborated, but I have nowhere else to go!
I am trying to fit some points to an inverse parabola of the form of $F(x)=1/(ax^2+bx+c)$.
My objective is to program a function in c++ that would take a set of 10-30 points and fit them to the inverse parabola.
I started trying to get an analytical expresion using Least squares, but I can't reach to get a result. I tried by hand (little crazy) and then I tried to solve analytically the expressions for a,b and c, but mupad doesn't give me a result (I am pretty new to Matlab's mupad so maybe i am not doing it correctly).
i don't know anymore how to approach the problem.
Can I get a anallytical expresion for this specific problem? I have also seen algorithms for general least squares fitting but I don't need a so complicated algorithm, I just need it for this equation.
If not, how would StackOverflow people approach the problem?
If needed I can post the equations, and the small Mupad code I've tried, but I think is unnecessary.
EDIT: some example
Im sorry the image is a little bit messy but it is the thing i need.
The data is in blue (this data is particularly noisy). I need to use only the data that is between the vertical lines (a bunch of data in the left and another one in the right).
The result of the fit is in de red line.
all this has been made with matlab, but I need to make it in c++.
I'll try to post some data…
Edit 2: I actually did the fitting in Matlab as follows (not the actual code):
create linear system Ax = b, with
A = [x² x 1]
x = [a; b; c]
b = 1/y;
It should work, shouldn it? I can solve then using Moore-Penrose pseudoinv calculated with SVD. Isn't it?
Best Answer
If you have some data $y$ and want to fit it to a curve of the form $1/(ax^2+bx+c)$, then what we can do is fit the data $1/y$ to the curve $ax^2+bx+c$. To do this in the least squares sense, we construct the following system of equations:
$$ \begin{pmatrix} x_1^2 & x_1 & 1 \\ x_2^2 & x_2 & 1 \\ \vdots & \vdots & \vdots \\ x_n^2 & x_n & 1\end{pmatrix} \begin{pmatrix} a \\ b \\ c\end{pmatrix} = \begin{pmatrix} 1/y_1 \\ 1/y_2 \\ \vdots \\ 1/y_n \end{pmatrix} $$
The matix on the left-hand side is known as a Vandermonde matrix. Denote this matrix $A$
Since in general $n \neq 3$, to solve this we simply solve $A^TA \mathbf{p} = A^T\mathbf{y}$, where $\mathbf{p} = (a,b,c)^T$ is the parameter vector.
Thus, $$\mathbf{p} = (A^TA)^{-1}A^T\mathbf{y},$$ which may be computed in many ways (note that $A^TA$ is symmetric, so the LU factorization should be easy to compute).
EDIT: Regarding largeness of $\mathbf{y}$ (which should not matter for the data in your picture).
Note that we are solving $A^TA\mathbf{p} = A^T\mathbf{y}$. Look at $A^T\mathbf{y}$:
$$A^T \mathbf{y} = \begin{pmatrix} \sum_{i=1}^n \frac{x_i^2}{y_1} \\ \sum_{i=1}^n \frac{x_i}{y_1} \\ \sum_{i=1}^n \frac{1}{y_1} \end{pmatrix}$$
If $y_i$ is so small that you are concerned about numerical stability, solve the following problem instead:
$$10^m \hat{a} x^2 + 10^m \hat{b} x + 10^m \hat{c} = 1/y$$
for some $m > 0$. Then, you condition your Vandermonde matrix in such a way that $A^T\mathbf{y}$ becomes:
$$A^T \mathbf{y} = \begin{pmatrix} \sum_{i=1}^n \frac{10^m x_i^2}{y_1} \\ \sum_{i=1}^n \frac{10^m x_i}{y_1} \\ \sum_{i=1}^n \frac{10^m}{y_1} \end{pmatrix}$$
and the smallness is gone. Compute $\hat{a},\hat{b},\hat{c}$ in the traditional sense, then $a = 10^m\hat{a}$, and so forth.