Message Boards Message Boards

GROUPS:

Minimizing a expression to determine a parameter

Posted 5 years ago
5383 Views
|
5 Replies
|
0 Total Likes
|

Hi, I have a complicated expression involving exponential terms. I want to find out the positive value of the only parameter for which the expression would be at minima. Use of "Minimize" operation results in the error message "Overflow occurred in computation" The expression is available in the attached file. thanks

Attachments:
POSTED BY: S G
Answer
5 Replies

Your expression involves high powers of the parameter. I think that's why you're getting overflow. Also, note that you're actually using NMinimize because your expression does not have infinite precision.

Posted 5 years ago

Frank is correct. One way to "fix" the problem is to replace D with d 10^(-20):

F = 1/3 (-11.620178705` + (
    1.2499999999999998` (1.9970632649872129` + E^(2.4137098720733767`*^20 (-6.0923968416000006`*^-24 - D))))/
    (1.996330429789108` + E^(3.0171373400917207`*^20 (-6.0923968416000006`*^-24 - D))) +
    (1.6666666666666665` (1.9970632649872129` + E^(2.4137098720733767`*^20 (-6.0923968416000006`*^-24 - D))))/
    (1.995110234550618` + E^(4.0228497867889607`*^20 (-6.0923968416000006`*^-24 - D))) +
    (2.4999999999999996` (1.9970632649872129` + E^(2.4137098720733767`*^20 (-6.0923968416000006`*^-24 - D))))/
    (1.9926743253237487` + E^(6.034274680183441`*^20 (-6.0923968416000006`*^-24 - D))))
F2 = Simplify[F /. D -> d 10^(-20)]
sol = Minimize[F2, d]
(* Solution in original units is *)
d 10^20 /. sol[[2]]

resulting in the following output:

 (E^(-2.41371 d) (-0.484169 E^(2.41371 d) + 0.0521019 E^(3.01714 d) + 
        0.0695118 E^(4.02285 d) - 0.864135 E^(5.43085 d) + 
        0.104395 E^(6.03427 d) - 0.829316 E^(6.43656 d) + 
        0.243228 E^(7.03999 d) - 0.759552 E^(8.44798 d) + 
        0.312996 E^(9.05141 d) - 1.45023 E^(9.4537 d) + 
        0.347816 E^(10.0571 d) - 1.3107 E^(11.4651 d) - 
        1.24106 E^(12.4708 d) + 0.904043 E^(13.0743 d) - 
        2.06531 E^(15.488 d)))/((0.499999 + 
        1. E^(3.01714 d)) (0.499998 + 1. E^(4.02285 d)) (0.499997 + 
        1. E^(6.03427 d)))
 {-3.87339, {d -> -269.7}}
 -2.697*10^22

A correction: While this gives you a result, there still is the issue of a loss of precision as when you plug in values surrounding this "optimal" value, you'll end up with the same (displayed and maybe "internally stored") value for the function. I'll have to think more about this.

Are you sure that your expression has a minimum? It seems to me that the infimum is attained as D tends to -Infinity. The structure of your expression is clearer if treated this way:

coeffs = Cases[F, _Real, All] // Union;
subst1 = Thread[coeffs -> Table[c[i], {i, Length[coeffs]}]];
c0 = (coeffs[[-4]]/12);
subst5 = Map[Reverse, 
   Cases[Rationalize[subst1[[1 ;; 9]]], Except[_Real -> _]]];
subst2 = Thread[
   Table[c[i], {i, 10, 13}] -> c[0]*Rationalize[coeffs[[-4 ;;]]/c0]];
subst3 = subst1 /. subst2;
F2 = Simplify[F /. subst3 /. subst5]
subst4 = c[0] (-D + c[2]) :> Log[x];
F3 = FullSimplify[F2 /. subst4]
F4 = F3 /. Map[Reverse, subst1]
Plot[F4, {x, -1.98^(1/15), 5}]

where x == Exp[c[0] (-D + c[2])] is a new variable.

Posted 5 years ago

Simple methods find one extremum, which is a maximum:

In[2]:= f = F /. D -> 1. 10^-24 d // Simplify;

In[3]:= df = D[f, d] // Simplify;

In[4]:= extremum = d /. (sol = Solve[df == 0, d, Reals][[1]])

During evaluation of In[4]:= Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result. >>

Out[4]= 3154.42

In[5]:= d2f = D[df, d] // Simplify;

In[6]:= (* extremum is a maximum *)
d2f /. sol

Out[6]= -2.05399*10^-8

In[7]:= (* with value *)
f /. sol

Out[7]= -1.88228

In[8]:= Plot[f, {d, extremum/10, 10 extremum}]

enter image description here

Posted 5 years ago

And

In[11]:= NMaximize[f, d]

Out[11]= {-1.88228, {d -> 3154.42}}

In[15]:= hl = Limit[f, d -> Infinity]

Out[15]= -2.06531

In[16]:= ll = Limit[f, d -> -Infinity]

Out[16]= -3.87339

In[23]:= Plot[f, {d, -50000, 50000}, 
 Prolog -> {Red, Line[{{-50000, ll}, {0, ll}}], 
   Line[{{0, hl}, {50000, hl}}]}, AxesLabel -> {"10^24 D", "F"}]

enter image description here

Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract