# Propagation of error through the solution for system of non-linear equations [closed]

I am currently working with a physical system governed by $n$ number of parameters. Thus the functional form of the behaviour of this system over a length-scale $z$ can be represented by a non-linear function $F(x_1,x_2,...x_n,z)$. Now, this physical system is divided into 2 sub-systems such that: $F(x_1,x_2,...x_n,z)=\begin{cases} F_1(x_1,x_2,...x_n,z) & \text{if$z< z_0$} \\ F_2(x_1,x_2,...x_n,z) & \text{if$z\geq z_0$} \end{cases}$

Out of these $n$ parameters, $n-2$ are known and were obtained through fitting various sets of data obtained through numerical simulations. Thus, there is an associated error $\delta x_i$ for parameters $x_i (i=1,2,...,n-2)$. A set of data also yields the numerical values of $F_1$ and $F_2$ at $z=z_0$. Let's call these values $f_1^{z_0}$ and $f_2^{z_0}$. Note that there is a discontinuity at $z_0$. These numerical values can be used to find the 2 remaining unknown parameters $x_{n-1}$ and $x_n$. For clarity, let's call them $y_1$ and $y_2$. Now the problem boils down to solving the following system of non-linear equations for $y_1$ and $y_2$:

$F_1(x_1,x_2,...x_{n-2},z_0,\mathbf{y_1},\mathbf{y_2})=f_1^{z_0}$

$F_2(x_1,x_2,...x_{n-2},z_0,\mathbf{y_1},\mathbf{y_2})=f_2^{z_0}$

This can be done easily by using fsolve in Matlab or FindRoot in Mathematica, and I obtain unique solutions for the 2 unknowns. However, I need to determine the error in the estimation of $y_1$ and $y_2$, as a result of error propagation from $x_1,x_2,...x_{n-2},f_1^{z_0},f_1^{z_0}$.

Does someone know of a systematic way to find $\delta y_1$ and $\delta y_2$?

So far I have tried the following 2 approaches:

1. A very crude method: Solve the pair of non-linear equations for combinations of 3 different values of each known parameter: $\{x_i-\delta x_i,x_i,x_i+\delta x_i\}$ and $\{f_1^{z_0}-\delta f_1^{z_0},f_1^{z_0},f_1^{z_0}+\delta f_1^{z_0}\}$, $\{f_2^{z_0}-\delta f_2^{z_0},f_2^{z_0},f_2^{z_0}+\delta f_2^{z_0}\}$. By doing so, I obtain $3^n$ sets of solutions for $y_1$ and $y_2$. Then I treat the $max\{y_1\}$ and $max\{y_2\}$ as my error bounds. I recently realised that the error propagation obtained in this way might be an overestimation.

2. Let $G_1=F_1(x_1,x_2,...x_{n-2},z_0,\mathbf{y_1},\mathbf{y_2})-f_1^{z_0}$ and $G_2=F_2(x_1,x_2,...x_{n-2},z_0,\mathbf{y_1},\mathbf{y_2})-f_2^{z_0}$. LEt the error for these functions be: $\delta G_1$ and $\delta G_2$. Now we simultaneously solve for $y_1$ and $y_2$ such that $G_1=0$ and $G_2=0$, $\implies$$\delta G_1=0$ and $\delta G_2=0$. Expanding the right-hand side of these equations using Taylor expansion up to first order by assuming that the errors are small enough, we get:

$\sum_{i=1}^{n-2}\frac{\partial G_1}{\partial x_i} \delta x_i+\frac{\partial G_1}{\partial f_1^{z_0}}\delta f_1^{z_0}+\frac{\partial G_1}{\partial y_1} \delta y_1+\frac{\partial G_1}{\partial y_2} \delta y_2=0$

$\sum_{i=1}^{n-2}\frac{\partial G_2}{\partial x_i} \delta x_i+\frac{\partial G_2}{\partial f_1^{z_0}}\delta f_1^{z_0}+\frac{\partial G_2}{\partial y_1} \delta y_1+\frac{\partial G_2}{\partial y_2} \delta y_2=0$

Then we can substitute the values of known parameters, their errors and solutions for $y_1$ and $y_2$ obtained by solving $G_1=0$ and $G_2=0$ in the above 2 equations we can find solutions for $\delta y_1$ and $\delta y_2$.

I have a feeling that this method assumes that the extraction of errors for all parameters is equivalent, which is not the case for $y_1$ and $y_2$, since they are evaluated using root-finding methods.

## closed as off-topic by sammy gerbil, Yashas, Kyle Kanos, Jon Custer, David HammenMar 11 '17 at 11:42

• This question does not appear to be about physics within the scope defined in the help center.
If this question can be reworded to fit the rules in the help center, please edit the question.

• Well, how did you find the error of the other parameters? Wasn't that a sort of root finding as well? – mikuszefski Mar 10 '17 at 14:57
• I found the errors on other parameters by fitting my data to different independent functions dependent on the individual parameters. Of course, non-linear least squares fitting is similar to root-finding, but there I had $m$ data-points (observations) and 1 unknown $(m>1)$ and the error is found in terms of the residual of the function. But now, we want to fit 2 data points $f_1^{z_0}$ and $f_2^{z_0}$ to 2 functions $F_1$ and $F_2$, in order to find $y_1$ and $y_2$. Thus, the residual will simply be the machine precision, or whatever accuracy I choose. – User1305 Mar 10 '17 at 16:17
• What I am looking for is an error estimate for $y_1$ and $y_2$ subject to uncertainties in all other parameters. – User1305 Mar 10 '17 at 16:19
• I'm voting to close this question as off-topic because this is a technical mathematical question with no physics content, and would be better on Mathematics SE or Cross-Validated SE. – sammy gerbil Mar 10 '17 at 17:38
• @sammygerbil While Cross-Validated SE might be the right place, one has too keep in mind: Error estimation, handling and propagation, concepts of errors, error handling, and statistics are of fundamental importance in experimental physics. – mikuszefski Mar 13 '17 at 6:54