Group Abstract Group Abstract

Message Boards Message Boards

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

How to choose starting point for FindMinimum

Cory Leahy
Posted 12 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
Sort By:
Cory Leahy
Cory Leahy
Posted 12 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
Cory Leahy
Cory Leahy
Posted 12 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.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard