# Mathematica gives complex coeffients when fitting polynomial

The problem is the following. I have two functions A[m,x] and B[m,x]. m is a variable, while x is a constant depending on other equations and parameters of the system. A is an integral which can only be evaluated numerically while B is a simple exponential function with the structure: B[m,x] = A[0,x]*Exp[m/C], where C is a constant (and so is A[0,x] for a given x).

Now, I fit a 20 degree polynomial to C = A[m,x]-B[m,x], in order to make Mathematica able to integrate C. This all appeared to work fine, until I suddenly discovered that for certain values of x, the polynomial fit has complex coefficients! This is very weird to me, because neither A nor B is complex. Could somebody please explain how the coefficients can become complex and how I can avoided it?

Thank you in advance!

EDIT: My code looks like basically this:

A[m_,x_] := NIntegrate[Exp[-(Sqrt[v^2 + 2*m] - V)^2/2], {v, 0, Sqrt[2*(x - m)]}]

B[m_,x_] := A[0,x] * Exp[m/T]

Tbcomb := Table[{mprim, N[A[mprim,x] - B[mprim,x]]}, {mprim, 0, x, 0.01}]
Combfit := Fit[Tbcomb,
{1, mprim, mprim^2, mprim^3, mprim^4, mprim^5, mprim^6, mprim^7, mprim^8,
mprim^9, mprim^(10), mprim^(11), mprim^(12), mprim^(13), mprim^(14),
mprim^(15), mprim^(16), mprim^(17), mprim^(18), mprim^(19),
mprim^(20)},
mprim
]


T and V are constants.

• Are you fitting or solving? Please post a self-contained functioning code. – Feyre Sep 7 '16 at 12:52
• If you're getting values like 10^-16 I then I would just Chop them as numerical artifacts. If they are sth like 0.5 I, then one can make assumptions in e.g. NonlinearModelFit: a\[Elem] Reals or so. – corey979 Sep 7 '16 at 12:55
• You can format inline code and code blocks by selecting the code and clicking the {} button above the edit window. The edit window help button ? is also useful for learning how to format your questions and answers. You may also find this meta Q&A helpful – Michael E2 Sep 7 '16 at 13:09
• @Feyre I am fitting. Code included! – Jhonny Sep 7 '16 at 13:10
• C is a reserved Mathematica symbol, so you should not use it as a variable. In fact, it is good practice to use only use variables names beginning with lower case letters in the Global context. – m_goldberg Sep 7 '16 at 13:10

## 1 Answer

Due to lack of a set of parameters V,T,x in the original question, maybe I'll solve some synthetic problem to illustrate the approach:

Let's generate some data from a 20th order polynomial with dispersion:

poly[y_] := 1 + Sum[y^i, {i, 1, 20}]
data = Table[{x, poly[x] + RandomReal[0.1]}, {x, -1, 1, 0.01}];
plot1 = ListPlot[data] The fit can be achieved with

f = Fit[data, Table[mprim^i, {i, 0, 20}], mprim]


1.03555 + 1.03503 mprim + 2.44307 mprim^2 - 2.50054 mprim^3 - 36.73 mprim^4 + 61.9126 mprim^5 + 412.885 mprim^6 - 468.637 mprim^7 - 2392.34 mprim^8 + 1986.07 mprim^9 + 8283.7 mprim^10 - 5017.65 mprim^11 - 17939. mprim^12 + 7804.28 mprim^13 + 24544.7 mprim^14 - 7313. mprim^15 - 20585.9 mprim^16 + 3789.38 mprim^17 + 9657.23 mprim^18 - 830.906 mprim^19 - 1937.02 mprim^20

We see that some of the coefficients are negative, but we want all of them to be positive.

Generate a polynomial func with coefficients coeff to be fitted:

var = Table[mprim^i, {i, 0, 20}];
coeff = Table[ToExpression["c" <> ToString[i]], {i, 0, 20}];
func = var.coeff;


and a list of conditions $c_i>0$:

cond = Table[coeff[[i]] > 0, {i, 1, Length@coeff}]


{c0 > 0, c1 > 0, c2 > 0, c3 > 0, c4 > 0, c5 > 0, c6 > 0, c7 > 0, c8 > 0, c9 > 0, c10 > 0, c11 > 0, c12 > 0, c13 > 0, c14 > 0, c15 > 0, c16 > 0, c17 > 0, c18 > 0, c19 > 0, c20 > 0}

Finally, one can use NonlinearModelFit:

nlm = NonlinearModelFit[data, {func, cond}, coeff, mprim];
fit = Normal[nlm]


1.04529 + 0.972808 mprim + 1.0053 mprim^2 + 1.20048 mprim^3 + 1.14997 mprim^4 + 0.0997689 mprim^5 + 0.688015 mprim^6 + 2.91753 mprim^7 + 0.20959 mprim^8 + 0.095588 mprim^9 + 2.86412 mprim^10 + 0.0702545 mprim^11 + 0.800427 mprim^12 + 0.114977 mprim^13 + 0.520907 mprim^14 + 2.81414 mprim^15 + 0.521147 mprim^16 + 1.69431 mprim^17 + 0.67438 mprim^18 + 0.0238567 mprim^19 + 1.56535 mprim^20

All of the coefficients are positive. To tackle the problem of constraining the coefficients of func to be real, in the definition of cond one ought to change ">0" into "\[Element] Reals", and in general it should work.