0
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ Are you fitting or solving? Please post a self-contained functioning code. $\endgroup$ – Feyre Sep 7 '16 at 12:52
  • 1
    $\begingroup$ 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. $\endgroup$ – corey979 Sep 7 '16 at 12:55
  • $\begingroup$ 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 $\endgroup$ – Michael E2 Sep 7 '16 at 13:09
  • $\begingroup$ @Feyre I am fitting. Code included! $\endgroup$ – Jhonny Sep 7 '16 at 13:10
  • 1
    $\begingroup$ 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. $\endgroup$ – m_goldberg Sep 7 '16 at 13:10
2
$\begingroup$

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]

enter image description here

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.

enter image description here

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.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.