Message Boards

WOLFRAM COMMUNITY

10845 Views

5 Replies

6 Total Likes

View groups...

Follow this post

Share this post:

GROUPS:

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 = (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:

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 12 years ago

Suggest you use FindRoot to solve your equations, rather than writing your own solver.

POSTED BY: Frank Kampas

Cory Leahy

Posted 12 years ago

FindRoot won't work for mutliple iterations at once. I had tried that originally

POSTED BY: Cory Leahy

Frank Kampas

Frank Kampas, Physicist at Large Consulting

Posted 12 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 BY: Frank Kampas

Shenghui Yang

Shenghui Yang, WOLFRAM

Posted 12 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"]]