Message Boards

WOLFRAM COMMUNITY

8954 Views

5 Replies

0 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

Material Sciences Physics

Minimizing a expression to determine a parameter

S G

Posted 9 years ago

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

5 Replies

Sort By:

David Keith

Posted 9 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.0539910^-8 In[7]:= (* with value *) f /. sol Out[7]= -1.88228 In[8]:= Plot[f, {d, extremum/10, 10 extremum}]

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 BY: David Keith

David Keith

Posted 9 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"}]

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

POSTED BY: David Keith

Gianluca Gorni

Gianluca Gorni, University of Udine

Posted 9 years ago

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.

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 BY: Gianluca Gorni

Jim Baldwin

Jim Baldwin, Retired

Posted 9 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.69710^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.

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.

POSTED BY: Jim Baldwin

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 9 years ago

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 BY: Frank Kampas

Reply to this discussion

Reply Preview

Attachments

Remove Add a file to this post

Follow this discussion

or Discard

Group Abstract

Feedback