4
$\begingroup$

I am trying to solve the FOC of this objective fcn oo with respect to x. I need to find optimal x with respect to k and l parameters that are seen in the objective.

However, although the integral and derivative look solvable, when I use Solve,NSolve, Reduce, Rationalize etc. I cannot solve this function and get x.

I receive an error:

Solve was unable to solve the system with inexact coefficients

So, do you have any idea how to fix this problem? P

This is the code:

pdf = PDF[NormalDistribution[a, b], y] /. {a -> 170497/10, b -> 121775/100}
 (2*Sqrt[2/Pi])/(E^((8*(-(170497/10) + y)^2)/23726641)*4871)
oo = x*(α + β*x) - Integrate[(x - y)*(α + β*y)*(1 - ℓ)*pdf, {y, 0, x}] + Integrate[(y - x)*(α + β*y)*(1 + \[ScriptK])*pdf, {y, x, 20845.3}]
x*(α + x*β) - (0.5000000000000001*(1. + \[ScriptK])*(E^(0.011497421822161848*x)*(-2.6341910136119725*^-40*α - 4.491216652478005*^-36*β) + 
     E^(3.37173728046882*^-7*x^2)*(-17010.989676334237*α + 0.9981723391216341*x*α - 2.913551090170966*^8*β + 17010.989676334237*x*β) + 
     E^(3.37173728046882*^-7*x^2)*((-17049.7 + 1.*x)*α + (-2.921751851525*^8 + 17049.7*x)*β)*
      Erf[9.900191736560279 - 0.0005806666238444241*x]))/E^(3.37173728046882*^-7*x^2) - 
  (1/(800*Sqrt[2*Pi]))*((-1 + ℓ)*((-97420*E^(58138454018/593166025)*(10*α + 170497*β) + 97420*E^((2*(170497 - 10*x)^2)/593166025)*
        (10*α + (170497 - 10*x)*β) + E^((4*(29069227009 - 1704970*x + 50*x^2))/593166025)*Sqrt[2*Pi]*
        ((6819880 - 400*x)*α + (116870074061 - 6819880*x)*β)*Erf[(170497*Sqrt[2])/24355])/
      E^((4*(29069227009 - 1704970*x + 50*x^2))/593166025) - Sqrt[2*Pi]*(40*(-170497 + 10*x)*α + (-116870074061 + 6819880*x)*β)*
      Erf[(Sqrt[2]*(-170497 + 10*x))/24355]))
deriv = D[oo, x]; 

Solve[deriv == 0, x]; 

During evaluation of In[20]:= PolynomialGCD::lrgexp:Exponent is out of bounds for function PolynomialGCD. >>

During evaluation of In[20]:= Solve::inex:Solve was unable to solve the system with inexact coefficients or the system obtained by direct rationalization of inexact numbers present in the system. Since many of the methods used by Solve require exact input, providing Solve with an exact version of the system may help. >>

FindRoot[deriv == 0, {x, 0}]

During evaluation of In[16]:= FindRoot::nlnum:The function value {0. +α-0.5 (<<1>>) (0. -1.77636*10^-46 (-17049.7 α-2.92175*10^8 β)+<<1>>+1. (0.998172 α+<<19>> β)+1. (1. α+17049.7 β))-0.000498678 (-1.+ℓ) (-4.45267*10^-46 (-6.81988*10^6 α-1.1687*10^11 β)+2.50663 (400. α+6.81988*10^6 β)+7.35018*10^-86 (3.41029*10^85 Plus[<<2>>]-3.59335*10^48 β-4.13142*10^45 Plus[<<2>>]-3.92096*10^83 Plus[<<2>>])+8.45081*10^-88 (5.15173*10^33 Plus[<<2>>]+3.41029*10^85 Plus[<<2>>]))} is not a list of numbers with dimensions {1} at {x} = {0.}. >>

FindRoot[deriv == 0, {x, 0}]
$\endgroup$
  • 1
    $\begingroup$ Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! $\endgroup$ – Michael E2 Aug 11 '16 at 15:30
  • $\begingroup$ People here generally like users to post code as Mathematica code instead of images or TeX, so they can copy-paste it. It makes it convenient for them and more likely you will get someone to help you. You may find this this meta Q&A helpful $\endgroup$ – Michael E2 Aug 11 '16 at 15:30
  • $\begingroup$ It looks like some code is missing (e.g. time, maybe more). Please clarify. $\endgroup$ – Michael E2 Aug 11 '16 at 15:32
  • $\begingroup$ time is the data that I use for fitting a distribution. I use time data and find the best fit according to FindDistribution, so that I get NormalDist with 17049.7 and 1217.75. The question that I asked is when I take the integrals and write the "oo" obj function and take partial derivative according to x, it looks fine,solvable. However, when I try to solve this partial derivative by Solve[derivative==0,x] I get that error message. $\endgroup$ – ecco Aug 11 '16 at 15:36
  • 1
    $\begingroup$ Then FindRoot[] is probably what you need instead of Solve[]. $\endgroup$ – Michael E2 Aug 11 '16 at 15:39
6
$\begingroup$

I do not think Mathematica can do exactly what you want it to do. There are two function you have tried, Solve and FindRoot, that come close; but neither is perfect.

Solve will take a given equation in x (say), along with some other parameters if desired, and attempt to find an exact solution for x in terms of whatever other parameters are present. It does this by a glorified version of how you or I would "solve for x" given an equation: it applies various manipulations and inverse functions to the equation in an attempt to isolate the desired variable. However, it cannot always do this, since not all equations can actually be solved in this way. As an example, the equation $x + \sin x = k$ has a single solution for $x$ for all values of $k$. But you can't write down an exact solution for $x$ in terms of "nice" functions like trig functions and their inverses, and so Mathematica can't use Solve to give you the value of $x$ as a function of $k$:

Solve[x + Sin[x] == k, x]
(* Solve::nsmet: This system cannot be solved with the methods available to Solve. >> *)

FindRoot, on the other hand, basically "looks at the graph of the function" and tries to figure out where the solution is. For example, if you want to solve the equation $x + \sin x = 1$, FindRoot will give you a solution:

FindRoot[x + Sin[x] == 1, {x, 0}]
(* {x -> 0.510973} *)

Internally, what Mathematica does here is the following:

  1. Start off at x = 0.
  2. Notice that the function $f(x) = x + \sin x$ is equal to zero at $x = 0$, and that it is increasing at this point.
  3. Reason that since we need to get to $f(x) = 1$, we should increase the value of $x$ and see what $f(x)$ is there.
  4. See what the function is doing at this new value of $x$. If $f(x)$ is too high, head in the direction that decreases $f(x)$. If it's too low, head in the direction that increase $f(x)$.
  5. Repeat step 4 until the value of $f(x)$ is sufficiently close to $f(x) = 1$, and return the result.

However, to do this procedure, you need to have a specific "target" value of $f(x)$ that you're aiming for. You can't just tell it to aim for a general value $k$, since step 4 above requires that you know whether $f(x)$ is "too high" or "too low" compared to the value you're aiming at. So you will not be able to solve the equation $f(x) = k$ for general $k$ in this way.

In your case, your function is a complicated expression in terms of four parameters, as well as a Gaussian function and the error function. I strongly suspect that these functions are not simple enough for the equation-manipulation algorithms of Solve (powerful though they are) to get an exact solution; and you cannot actually use FindRoot without having specific values of the coefficients in mind.

You can, if you want, use FindRoot to solve for specific values of your parameters; for example:

FindRoot[(deriv /. {\[Alpha] -> 0, \[Beta] -> 3, \[ScriptK] -> 0, \[ScriptL] -> 0.3}) == 0, {x, 0}]
(* {x -> 8515.17} *)

The /. construction temporarily sets the parameters of your problem to the specific values given; and then FindRoot can actually figure out where the solution for the equation deriv == 0 is for these specific parameters.

EDIT: If you want to automate this to create a table of solutions for various values of $k$ and $l$ while keeping $\alpha$ and $\beta$ fixed, you can use Table to automate this somewhat:

x /. Table[ FindRoot[(deriv /. {\[Alpha] -> 0, \[Beta] -> 3, \[ScriptK] -> k0, 
    \[ScriptL] -> l0}) == 0, {x, 0}], {k0, 0, 1, 0.2}, {l0, 0, 1, 0.2}]

This returns a 6x6 table of solutions for $x$ for the fixed values of $\alpha = 0$ and $\beta = 3$, and the varying values $(k,l) = (0,0), (0, 0.2), (0, 0.4), ..., (0.2, 0), (0.2, 0.2), ...$ A more "finely spaced" table could be created by changing the 0.2 in the above command to a smaller number (it corresponds to the increment between successive $k$ and $l$ values.

| improve this answer | |
$\endgroup$
  • $\begingroup$ Thank you for your effort, it is really helpful to understand the working principle of these two functions. If I take alpha and beta parameters as fixed, is there a shorter way (like a for loop) to find all the x values for the combinations of k and l parameters which can take any value between 0 and 1? $\endgroup$ – ecco Aug 11 '16 at 17:42
  • $\begingroup$ Not all possible values, but you could make a table of solutions. I'll edit my answer to show how this could be done. $\endgroup$ – Michael Seifert Aug 11 '16 at 17:54
  • $\begingroup$ Related to this table issue,how can I make a graph of the calculated values in the table? I am trying to construct a 3D graph with two parameters and a corresponding solution. Thank you $\endgroup$ – user42251 Aug 11 '16 at 18:24
  • $\begingroup$ @bbusra: Look up ListPlot3D. $\endgroup$ – Michael Seifert Aug 11 '16 at 18:58

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.