Message Boards

GROUPS:

Cory Leahy

How to choose starting point for FindMinimum

Cory Leahy

Posted 10 years ago

9814 Views

5 Replies

6 Total Likes

Follow this post

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 = (RT/v^2) (1 - c/(vT^3)) (v + B (1 - b/v)) - (A (1 - a/v))/v^2 - P; FXX = (1/(T^2v^5)) ((-aT^2vA (1 - a/v)) + 2T^2v^2A (1 - a/v) + R (bB (3T^3v - 4 c)) + (v((3Bc) - (2BT^3 v) + (2cv + 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

Posted 10 years ago

Wow, thanks for your help. This was really informative!

POSTED BY: Cory Leahy

Shenghui Yang

Shenghui Yang, WOLFRAM

Posted 10 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 control <Omit the definition of constants> FX[v_, P_] := (RT/v^2) (1 - c/(vT^3)) (v + B (1 - b/v)) - (A (1 - a/v))/v^2 - P; FXX[v_]:=(1/(T^2v^5)) ((-aT^2vA (1-a/v))+2T^2v^2A (1-a/v)+R (bB (3T^3v-4 c))+(v((3Bc)-(2BT^3 v)+(2cv+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

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 10 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}}