Message Boards Message Boards

GROUPS:

How to choose starting point for FindMinimum

Posted 10 years ago
9814 Views
|
5 Replies
|
6 Total Likes
|
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 10 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 10 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