3
$\begingroup$

I have the fallowing code (the functions within that code are replaced by simpler ones, but the code structure is identical): $$$$

g[x_, y_, z_] := x^2 + y^2 + z^2;
f[x_, y_, z_] := (buffer = Minimize[g[x, y, z], x];
   xmin = x /. buffer[[2]];
   g[xmin, y, z]);

Plot3D[f[x, y, z], {y, 0, 4}, {z, 0, 4}, PlotPoints -> 2]


It takes my computer about 7 seconds to process that code. With my actual functions it takes several minutes. This is only one step in a longer chain of steps my program needs to process.

Is there a way to speed that step up? Maybe something faster than using the function Minimize.

$$$$ $$$$ If i modify my code like this:

 g[x_, y_, z_] := x^2 + y^2 + z^2;
 f[x_, y_, z_] := (buffer = Minimize[g[x, y, z], x];
       xmin = x /. buffer[[2]];
       g[xmin, y, z]);    
Plot3D[f[x, y, z], {y, 0, 4}, {z, 0, 4}, PlotPoints -> 2,
             RegionFunction -> Function[{x, y, z}, x + y < 3]]

$$$$

The result is what i want, but it actually takes more time, even though in the resulting plot are less data points. I assume mathematica plots first my function and then it cuts with RegionFunction in my plot. Is there a way to make Mathematica just plot the resulting area and not plot the area that is later cut out anyway?

I've tried it like this:

 g[x_, y_, z_] := x^2 + y^2 + z^2;
 f[x_, y_, z_] := (buffer = Minimize[g[x, y, z], x];
       xmin = x /. buffer[[2]];
       g[xmin, y, z]);    
Plot3D[f[x, y, z], {y, 0, 4}, {z, 0, 4 - y}, PlotPoints -> 3]

$$$$

But the resulting plot doesn't exactly look like the one resulting using RegionFunction, it looks abit more 'bumpy'. So i've increased the PlotPoints from 2 to 3, then the resulting plot looks pretty much like the one with RegionFunction. But the Problem is, now plotting takes longer then with RegionFunction ;)

$$$$ Thanks for your help.

$$$$ $$$$ Your answeres (using the function Evaluate[]) helped me solve the problem for the above described situation. However it does not work with my actual functions, so i have to write them down here after all.

bindeaa[e0_] := e0;
c[e0_] := 40 e0;
bindeab[e0_] := 1.27 e0; 

alpha[gamma_, p_, epsilon_] := p (gamma  - epsilon);
beta[b_, epsilon_, n2_] := b epsilon^2 n2^(3/2);  


eisl[b_, e0_, g_, gamma_, p_, epsilon_, n2_, x_] := 
            g c[e0] (epsilon)^2 - bindeaa[e0] + 
            e0 ((-2 Log[E^(1/2) x])/x^2 + alpha[gamma, p, epsilon]/x + 
            beta[b, epsilon, n2]/x^(3/2));


eislminx[b_, e0_, g_, gamma_, p_, epsilon_, n2_]:= 
(minxzwischen = FindMinimum[{eisl[b, e0, g, gamma, p, epsilon, n2, x], 0 < x},x];
              minx = x /. minxzwischen[[2]];          
              eisl[b, e0, g, gamma, p, epsilon, n2, minx]);

 Plot[eislminx[10, 1, 0.7, 0.3, 4.9, 0.05, n2], {n2, 0, 10}, 
 PlotPoints -> 2, Evaluated -> True]

Mathematica gives the following error messages:

FindMinimum::nrgnum: The gradient is not a vector of real numbers at {x} = {1.}. >> FindMinimum::nrgnum: The gradient is not a vector of real numbers at {x} = {1.}. >> FindMinimum::nrgnum: The gradient is not a vector of real numbers at {x} = {1.}. >> General::stop: Further output of FindMinimum::nrgnum will be suppressed during this calculation. >> ReplaceAll::reps: {x} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >> ReplaceAll::reps: {x} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing. >>

What is the problem here?

Improve this question
$\endgroup$
3
  • $\begingroup$ Your code is not complete so your example cannot be evaluated. Please correct that, and also format your code using the StackExchange editor. $\endgroup$
    – Mr.Wizard
    Jun 27, 2014 at 7:43
  • $\begingroup$ I've edited my question and i think it fit's the rules now. $\endgroup$
    – 11drsnuggles11
    Jun 28, 2014 at 9:09
  • $\begingroup$ What is missing in my code? i can process the code as it is and get the plot i expect to get. $\endgroup$
    – 11drsnuggles11
    Jun 28, 2014 at 9:41

2 Answers 2

Reset to default
4
$\begingroup$

We can find the minimum-producing argument x as a function of n2 using NDSolve. We need only find an initial minimum and a differential equation for the trace of the minima. In fact, we can integrate the minimum value of eisl and return a single function that yields the minimum for each n2.

Given that the parameters b,e0,g,gamma,p,epsilon`, are (constant) numbers, the differential equation for the minima may be obtained from

Dt[D[eisl[b, e0, g, gamma, p, epsilon, n2, x], x]] == 0 /. x -> x[n2] /. Dt[n2] -> 1

If the position n2, x is at a minimum, then D[eisl[b, e0, g, gamma, p, epsilon, n2, x] will be zero. Then as n2 increases, we want n2, x to remain at a minimum. Therefore, partial derivative should not change value but remain at zero. Therefore, the differential

Dt[D[eisl[b, e0, g, gamma, p, epsilon, n2, x], x]]

should be zero. Hence the differential equation above. This will yield x as a function x[n2] of n2.

To integrate the minimum value of eisl, we have NDSolve integrate

y'[n2] == D[eisl[b, e0, g, gamma, p, epsilon, n2, x[n2]], n2]

which can be accomplished at the same time as NDSolve integrates x[n2].

There is one difficulty. For the parameter values in the OP's example, there is no minimum when n2 is greater that 4.97 or so. For values of n2 greater than this, the function eisl decreases to a limit -0.93 as x -> Infinity. I will assume that the limit should be the value of the function. To handle that exception, I will combine the solution from NDSolve with the limit using Piecewise.

Aside: The real value may be found with the following code.

FindRoot[{D[eisl[10, 1, 0.7, 0.3, 4.9, 0.05, n2, x], x] == 0, 
  D[eisl[10, 1, 0.7, 0.3, 4.9, 0.05, n2, x], {x, 2}] == 0}, {x, 
  1}, {n2, 4.9}]
(*
  {x -> 2.97292, n2 -> 4.96976}
*)

To get a little closer to the boundary value, I used the option Method -> "StiffnessSwitching".

Code

Below is a rewritten eislminx function. It returns a pure Function, so you have to pass it an argument (an input value for n2).

eislminx[
    b_?NumericQ, e0_?NumericQ, g_?NumericQ, gamma_?NumericQ, p_?NumericQ, epsilon_?NumericQ] := 
 eislminx[b, e0, g, gamma, p, epsilon] = (* memoize solution *)
  Module[{y, x, x0, sol, n2},
   x0 = First@
     FindArgMin[{eisl[b, e0, g, gamma, p, epsilon, 0, x], 0 < x}, x];
   sol = NDSolveValue[{
      y'[n2] == D[eisl[b, e0, g, gamma, p, epsilon, n2, x[n2]], n2],
      Dt[D[eisl[b, e0, g, gamma, p, epsilon, n2, x], x]] == 0 /. x -> x[n2] /. Dt[n2] -> 1,
      x[0] == x0, y[0] == eisl[b, e0, g, gamma, p, epsilon, 0, x0]},
     y, {n2, 0, 10}, Method -> "StiffnessSwitching"];
   With[{
     sol1 = sol,
     lim = Limit[eisl[b, e0, g, gamma, p, epsilon, n2, x],
       x -> Infinity, Assumptions -> n2 > sol["Domain"][[1, 2]]],
     n1 = sol["Domain"][[1, 2]]}, 
    Function[{n}, Piecewise[{{sol1[n], 0 <= n <= n1}}, lim]]
    ]
   ]

Example plot: The first time eislminx is called there will be a warning message indicating that NDSolve ran into the value of n2 where there ceases to be a minimum. Ignore it, or use Quiet.

Plot[eislminx[10, 1, 0.7, 0.3, 4.9, 0.05][n2], {n2, 0, 10}] // AbsoluteTiming

NDSolveValue::ndsz: At n2$130831 == 4.969743353334403`, step size is effectively zero; singularity or stiff system suspected. >>

Mathematica graphics

The second time you run Plot, it will be much faster, because the solution of the differential equation has been saved (for the particular values of the parameters). See What does the construct f[x_] := f[x] = ... mean?.

Plot[eislminx[10, 1, 0.7, 0.3, 4.9, 0.05][n2], {n2, 0, 10}]; // AbsoluteTiming
(*
  {0.022850, Null}
*)
Improve this answer
$\endgroup$
1
  • $\begingroup$ @11drsnuggles11 You're welcome. :) $\endgroup$
    – Michael E2
    Jul 2, 2014 at 12:25
2
$\begingroup$

Regarding the updated version of your question, the problem is that there is still some attempt being made inside Plot to evaluate the function symbolically, leading to the error messages from FindMinimum requiring a numerical starting point, not a symbolic one.

This kind of thing happens a lot, and the most common way to fix it is to define the function only for numerical values of the plot variable. This is done by adding ?NumericQ to the pattern for the last variable in eislminx:

ClearAll[eislminx];
eislminx[b_, e0_, g_, gamma_, p_, epsilon_, n2_?NumericQ] := 
 Module[{minxzwischen, minx, x},
  (minxzwischen = 
    FindMinimum[{eisl[b, e0, g, gamma, p, epsilon, n2, x], 0 < x}, 
     x];
   minx = x /. minxzwischen[[2]];
   eisl[b, e0, g, gamma, p, epsilon, n2, minx])]

Plot[eislminx[10, 1, 0.7, 0.3, 4.9, 0.05, n2], {n2, 0, 10}, 
 MaxRecursion -> 1, PlotRange -> {-.93, -.82}]

plot

This makes sure that that eislminx will not be entered unless the variable n2 has a numerical value.

In Plot, I also replaced the choice of PlotPoints (which only provides a starting value for the mesh) by a setting for MaxRecursion. Alternatively, when I want more control over the points at which to evaluate the plot function, it's sometimes better to do this:

ListLinePlot[
 Table[eislminx[10, 1, 0.7, 0.3, 4.9, 0.05, n2], {n2, 0, 10, .1}], 
 PlotRange -> {-.93, -.82}, DataRange -> {0, 10}]

which produces a plot that is indistinguishable from the above, but is based on a Table of rigidly defined values.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thanks, Jens. But your code isn't speeding up things compared to my last posted code without the function "Evaluate", your's is actually slower. Maybe i expressed myself in an unclear way, i tried to make that code work with the function Evaluate, because this function speeded up the process of the initial, simplified code sample, so i hoped it would do the same with the code at the end of my question. Without evaluate that code at the end of my question works fine, but is slow...but i don't know how to make it work with Evaluate $\endgroup$
    – 11drsnuggles11
    Jun 29, 2014 at 8:35
  • $\begingroup$ I think my code is slower than your (without Evaluate) because I replaced PlotPoints by MaxRecursion. But I can't think of a way to make it faster (even though you could add Evaluate to my solution, the speed doesn't change). I thought you wanted to get rid of the error messages, and using ?NumericQ is the most speed-friendly way I know that does this. $\endgroup$
    – Jens
    Jun 29, 2014 at 17:11
  • $\begingroup$ On my computer, the following is faster in the Plot command: PlotPoints -> 8, MaxRecursion -> 2 $\endgroup$
    – Jens
    Jun 29, 2014 at 17:22
  • $\begingroup$ Yes, i wanted to get rid of the error messages, and thanks again for that, Jens. $\endgroup$
    – 11drsnuggles11
    Jun 29, 2014 at 18:47

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.