Message Boards Message Boards

1
|
10330 Views
|
5 Replies
|
6 Total Likes
View groups...
Share
Share this post:

How to choose starting point for FindMinimum

Posted 11 years ago
My goal is to use the Newton Raphson root finding method with the Ideal Gas Law, Van der Waal equation and Beattie-Bridgeman equation to solve for the Specific Volume over a process of constant Temperature but variant Pressure.

I was able to get the Van der Waal equation, assuming the code should work for the other equations, applied them.  However, as I manually worked out the Root Finding Method for the first few iterations for all equation, I know what is expected I'm not getting the same results for the Ideal Gas Law nor the Beattie-Bridgeman Equation. In fact I'm getting very odd results with some errors and I suspect it must have to do with what starting number I choose for the FindMinimum function.

This is my code Beattie-Bridgeman:
 R = 8.314;
 A = 507.28;
 B = 0.10476;
 a = 0.07132;
 b = 0.07235;
 c = 660000;
 T = 333;
 iter = 60;
 
FX = (R*T/v^2) (1 - c/(v*T^3)) (v + B (1 - b/v)) - (A (1 - a/v))/v^2 - P;

FXX = (1/(T^2*v^5)) ((-a*T^2*v*A (1 - a/v)) + 2*T^2*v^2*A (1 - a/v) + R (b*B (3*T^3*v - 4 c)) + (v*((3*B*c) - (2*B*T^3 v) + (2*c*v + T^3)*(-(v^3)))));

FXXX = v - FX/FXX;

nn = Table[FXXX, {P, 400, 6000, iter}];

Ps = Table[P, {P, 400, 6000, iter}];

vs = Table[
]  v /. Last@FindMinimum[nn[[i]], {v, 7 i^(-.5) - 0.5}], {i, Length@nn}]

ListPlot[Thread[{Ps, vs}], Joined -> True, PlotRange -> {0, 7}]


Output:
 {6.83167, 8.25229, 9.92447, 11.6632, 13.3806, 15.0489, 16.6622, \
 ]18.2218, 19.7311, 21.1941, 1.61453, 1.61453, 1.61453, 1.61453, \
 1.61453, 1.61453, 1.61453, 1.14992, 1.61453, 1.61453, 35.0641, \
 1.33489, 1.18791, 1.08221, 0.996586, 0.923465, 0.859006, 0.800873, \
 0.747449, 0.697476, 0.649823, 0.603252, 0.555946, 0.503696, -1.76313, \
 -1.76313, 0.306443, -1.76313, -1.76313, -1.76313, 0.196067, -1.76313, \
 -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, \
 -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, \
 -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, \
-1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, \
-1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, \
1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, -1.76313, \
-1.76313, -1.76313, -1.76313, 0.246203, -1.76313, -1.76313, -1.76313, \
-1.76313, -1.76313, -1.76313}

and

Errors:
FindMinimum::lstol: The line search decreased the step size to within the tolerance specified by AccuracyGoal and PrecisionGoal but was unable to find a sufficient decrease in the function. You may need more than MachinePrecision digits of working precision to meet these tolerances. >>

General::stop: Further output of FindMinimum::nrnum will be suppressed during this calculation. >>

Power::infy: Infinite expression 1/0. encountered. >>

General::stop: Further output of FindMinimum::nrnum will be suppressed during this calculation. >>


Any help would be greatly appreciated.  Thanks!
POSTED BY: Cory Leahy
5 Replies
Posted 11 years ago
Wow, thanks for your help.  This was really informative!
POSTED BY: Cory Leahy
@Cory, 

To solve very complicated and non-linear equations, I usually visualize them before I use numerical solvers or optimization function. In your case, the same rule applies. Also, if you use function instead of algebraic expression, you can have better control
<Omit the definition of constants>
FX[v_, P_] := (R*T/v^2) (1 - c/(v*T^3)) (v +
B (1 - b/v)) - (A (1 - a/v))/v^2 - P;
FXX[v_]:=(1/(T^2*v^5)) ((-a*T^2*v*A (1-a/v))+2*T^2*v^2*A (1-a/v)+R (b*B (3*T^3*v-4 c))+(v*((3*B*c)-(2*B*T^3 v)+(2*c*v+T^3)*(-(v^3)))));
FXXX[v_,p_]:=v-FX[v,p]/FXX[v];
Now you can get rid of nn and use FXXX function directly (with P = 400)

To check the local minimum for a wide range of P, Manipulate function is rather handy:
Manipulate[Plot[FXXX[v,P],{v,-10,10}],{P,400,6000,iter}]

It seems that there are singularities over the (-Inf, Inf) domain. If you are just interested in the minima over x<-2  (approx. ) , you just need: 
(res=FindMinimum[{FXXX[v,#],v<-2},{v,-5}]&/@Range[400,6000,iter])//TableView

If I extract the result from the above (vs = all value of v from the table above, ps = "Ps" in your code and mins= the first colum in the table above), I have a manifold of minima in terms of vs and ps:
vs=v/.res[[All,2]];
ps=Range[400,6000,iter];
mins=res[[All,1]];
ListPointPlot3D[
      Transpose[{vs,ps/1000,Log[10,mins]}],
      Filling->Bottom,AxesLabel->{"v","ps/1000\n","Log10[Mins]    "},
      LabelStyle-> Directive[16,Bold,FontFamily->"Helvetica"],
      BoxRatios->{2,1,1}
]

 

Adding some extra Mathematica graphical function can make this even better. We can easily put the discrete plot and a smooth surface of FXXX together with Show function. The options in plt1 and plt2 (such as a->b)  are just for artistic purpose and you can live without them if you just want to run the simple code. 
minimaPoints = Transpose[{vs, ps/1000, Log[10, mins]}]
plt1=ListPointPlot3D[minimaPoints,FillingStyle-> Orange,PlotStyle-> {{White,PointSize[0.015]}},Filling->Bottom];
plt2=Plot3D[Log[10,FXXX[v,1000*p]],{v,-8,-2.1},{p,0.4,6},PlotStyle->Directive[Opacity[0.7],Specularity[1,20]], ColorFunction->"SunsetColors" , Mesh->None];
Show[plt1, plt2, BoxRatios -> {0.8, 1, 0.4}, Boxed -> False, PlotRange -> {{-8, -2.1}, {0.4, 6}, Automatic},  AxesLabel -> {"v", "ps/1000", "Log10[Mins] "},  LabelStyle -> Directive[16, Bold, FontFamily -> "Helvetica"]]

POSTED BY: Shenghui Yang
Not sure I know what you mean by "multiple iterations at once".  You can do something like this:
In[1]:= FindRoot[ Sin[x ]== 0, {x,#}] &  /@  Range[5]

Out[1]= {{x->0.}, {x->3.14159}, {x->3.14159}, {x->3.14159}, {x->9.42478}}
POSTED BY: Frank Kampas
Posted 11 years ago
FindRoot won't work for mutliple iterations at once.  I had tried that originally
POSTED BY: Cory Leahy
Suggest you use FindRoot to solve your equations, rather than writing your own solver.
POSTED BY: Frank Kampas
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