1
|
10396 Views
|
5 Replies
|
6 Total Likes
View groups...
Share
GROUPS:

# 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}andErrors: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!
5 Replies
Sort By:
Posted 11 years ago
 @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 controlFX[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])//TableViewIf 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 11 years ago
 Suggest you use FindRoot to solve your equations, rather than writing your own solver.
Posted 11 years ago
 FindRoot won't work for mutliple iterations at once.  I had tried that originally
Posted 11 years ago
 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 11 years ago
 Wow, thanks for your help.  This was really informative!