Why Would fplot(f) and fplot(vpa(f)) Show Different Results?

I'm seeing different results between fplot for a symfun object, the fplot of the vpa form of that object, and the regular old plot after evaluating and casting it to a double. Is there something about fplot that I'm missing or is this an issue peculiar to this symfun? Other than that erfi at the end of the function, I don't see anything particularly striking about this function. I'm pretty sure that the vpa form and the cast to double are correct..

>> whos fu1
  Name      Size            Bytes  Class     Attributes

  fu1       1x1                 8  symfun              

>> fplot(fu1,[0 15])
>> hold on
>> plot(svec,double(fu1(svec)),'r')
>> fplot(vpa(fu1),[0 15],'x')
>> fu1
 
fu1(y1) =
 
piecewise(y1 < 0, 0, 0 <= y1, -(2535301200456458802993406410752*exp(-(2*y1)/5)*exp(-(3*y1)/10)*exp(-(y1/10 + 1)^2/2)*exp(-(y1/5 + 1)^2/8)*abs(y1)^3)/(31859534503965572279823959492121*((1053932631846001350838118399344640000*pi^(1/2)*exp(279/16))/31859534503965572279823959492121 - (921805937454099571820119166106277959316698824704*exp(-6250000000850000000001/10000000000000000000000))/121534479156362809294982755630954742431640625 + (pi^(1/2)*exp(279/16)*erfi(425000000001i/100000000000)*1053932631846001350838118399344640000i)/31859534503965572279823959492121)))

>> vpa(fu1)
 
ans(y1) =
 
piecewise(y1 < 0, 0.0, 0 <= y1, 0.20729535137578417577792607729466*y1^3*exp(-0.125*(0.2*y1 + 1.0)^2)*exp(-0.5*(0.1*y1 + 1.0)^2)*exp(-0.3*y1)*exp(-0.4*y1))

[Walter Roberson]

It is common to see differences between the two when roots are involved, as high degree roots are a real balancing act. But this particular situation does not have roots so I would need further investigation.

[Paul]

Interesting. I have another example as well that doesn't involve roots. But in that other example the fplot is clearly wrong (whereas here it looks like it could be correct) and fplot(vpa()) is also clearly wrong.

Also, FWIW the red and blue curves differ by a factor of ~1.53

Answer 1

[Walter Roberson]

fplot problems. Look at this:

format long g
syms y1
fu1 = str2sym("piecewise(y1 < 0, 0, 0 <= y1, -(2535301200456458802993406410752*exp(-(2*y1)/5)*exp(-(3*y1)/10)*exp(-(y1/10 + 1)^2/2)*exp(-(y1/5 + 1)^2/8)*abs(y1)^3)/(31859534503965572279823959492121*((1053932631846001350838118399344640000*pi^(1/2)*exp(279/16))/31859534503965572279823959492121 - (921805937454099571820119166106277959316698824704*exp(-6250000000850000000001/10000000000000000000000))/121534479156362809294982755630954742431640625 + (pi^(1/2)*exp(279/16)*erfi(425000000001i/100000000000)*1053932631846001350838118399344640000i)/31859534503965572279823959492121)))")

fu1 = 

G = simplify(fu1,'steps', 50)

G = 

%G is much more compact than fu1, but is it equal?

fplot(fu1, [0 15])
title('fu1')

%... its is empty except for a possible discontinuity ??

%... was it complex valued maybe?
fplot(real(fu1), [0 15])
title('real(fu1)')

%nope, that doesn't show anything

%now let us look at the compact version:

fplot(G, [0 15])
title('G')

%... it is all negative???

%did we get some kind of complication from it being complex?
fplot(real(G), [0 15])
title('real(G)')

%nope, not visibly.

%so take a sample value, see if we get complex or positive or negative
double(subs(fu1,y1,1))

ans = 
        0.0469528822979662

vpa(subs(fu1,y1,1))

ans = 

0.046952882297966167119634566815625

double(subs(vpa(fu1),y1,1))

ans = 
        0.0469528822979662

%those are all the same and all positive, so details of the processing
%order should now matter much

%how about for the compacted
double(subs(G,y1,1))

ans = 
        0.0469528822979662

vpa(subs(G,y1,1))

ans = 

0.046952882297966167119634566875553

double(subs(vpa(G),y1,1))

ans = 
        0.0469528822979662

%those are all the same and all positive, and the same as fu1.
%so at least at this test point, the compacted version is the
%same as the longer version

%how about numeric plotting?

X = linspace(0,15);
Y = double(subs(fu1, y1, X));
plot(X, Y)
title('fu1 numeric')

Y2 = double(subs(G, y1, X));
plot(X, Y2)
title('G numeric')

%those are the same, and positive, and the shape you are expecting

So fplot is not acting reasonably here.

[Walter Roberson]

I filed a bug report pointing to here.

[Paul]

Thanks for this repsonse.

I had tried simplify too but not with as many steps. It looks like your plots "fu1 numeric" and "G numeric" match the red curve in my original question. But why do you think that none of your other plots matched the blue curve in my original question? I suspect it has something to do with the fact that str2sym doesn't do a literal copy of the argument to the LHS. It looks like it does some simplifying on its own. Look at fu1 as compared to the input to str2sym.

And then I don't know what to think about plots of G and real(G).

Because the bug report is pointing to this thread I may come back later and include a second example, though as you'll see it's related to the first.

I trust that when you get a response from Tech Support on the bug report that you'll let us know here what they say.

Thanks again.

[Paul]

Walter,

As you can see the function fu1 is a ratio and was, in fact, formed as the explicit ratio of two other functions, e.g.,

fu1 = g2/g1;

As a general rule, is it better to simplify or vpa before combining functions algebraically? In this particular case I checked and

fu1 = vpa(g2)/vpa(g1);

Seemed to woek better in this case, just wondering if that's best practice in general.

[Walter Roberson]

I only recommend vpa() as a final step, not as an intermediate step. Even as a final step it can be risky.

On the other hand, it can be useful to simplify a numerator and denominator separately.

[Paul]

@Walter Roberson,

Did you get a response from Tech Support?

[Walter Roberson]

It was a Thanks for the bug report kind of answer.

[Christopher Creutzig]

There is a problem with some intermediate calculation being done at low precision, compare:

digits(5)
X = linspace(0,15);
Y2 = double(vpa(subs(G, y1, X)));
plot(X, Y2)
title('G numeric @ digits==5')

Since that problem affects a constant in G, fplot(vpa(G),[0 15]) seems like a reasonable workaround.

[Walter Roberson]

I have had too many cases that mathematically exploded into roots of higher order polynomials involving large coefficients. vpa() of such systems are often mathematical disasters as the delicate balances between components are not kept.

I have also had too many cases involving bessel functions that were ruined by premature vpa() :(

[Paul]

@Christopher Creutzig

Can you expand on your comment? Can you say where the problem lies? i ran your code and got an incorrect plot that looked very close (if not identiical) to Walter's plots titled G and real(G), which I suspect was the point of your comment.

My bigger concern is with the blue curve in my original question, which I'm quite certain is wong but is not obviously wrong. I mean, if fplot shows nothing it's a clue something might be wrong. But that blue curve is at least plausible.

Also, do you know what str2sym does wrt to simplification? Apparently it doesn't do a lteral copy, which I suspect has something to do with that blue curve not being reprodcuible here.

[Christopher Creutzig]

@Paul

The problem is in the code fplot runs to evaluate symbolic expressions. It will be fixed in a future release.

I agree that the plot not being obviously wrong is a problem. The background as to what happens here is that the function essentially is complicatedConstant*f(x), and vpa(complicatedConstant) with low precision returns an incorrect value. That incorrect value then scales all function evaluations.

As you may know, sufficiently complicated expressions cannot be converted to floating point numbers correctly, no matter what - same class of "cannot be done in general" as the halting problem. Which means we (just like everyone else) can only mitigate the issue, not solve it once and for all.

str2sym does the same simplifications as if you typed in the symbolic expression, declaring symbolic variables, constants, etc:

str2sym("1 + 1") == sym(1) + sym(1)     % aka sym(2)
str2sym("x - x") == sym('x') - sym('x') % aka sym(0)
str2sym("sin(x)^2 + cos(x)^2") == sin(sym("x"))^2 + cos(sym("x"))^2 % does not auto-simplify to 1

[Christopher Creutzig]

@Walter Roberson Agreed, I do not recommend to run vpa() on intermediate results in general. The examples you gave are some of the reasons for that.

[Walter Roberson]

There are cases for symbolic work where what gets displayed involves different names than exist internally. Typing those expressions in as text (like I did) would then not get back the same expression.

The situation is sort of like writing, "John goes 5 feet left and John goes 3 feet right" and concluding that John ends up 2 feet left, when internally it was "John (Lennon)" in one case and "John (Prine)" in the other case.

[Paul]

@Christopher Creutzig

Thanks for the response.

Is the issue with fplot a formal bug? if so, can you post the link to the bug report?

I'm confused about this vpa(complicatedConstant) explanation. This command:

fplot(fu1,[-1 10])

shows (2019a) a plot that is zero from -1 to 0, a vertical dashed line indicating a singularity at 0, and then nothing from 0 to 10. So even if fplot is using vpa with low precision under its hood, the problem doesn't seem to be just a simple scaling across the entire plot domain.

Also, conider these results:

clear
fu1 = str2sym("piecewise(y1 < 0, 0, 0 <= y1, -(2535301200456458802993406410752*exp(-(2*y1)/5)*exp(-(3*y1)/10)*exp(-(y1/10 + 1)^2/2)*exp(-(y1/5 + 1)^2/8)*abs(y1)^3)/(31859534503965572279823959492121*((1053932631846001350838118399344640000*pi^(1/2)*exp(279/16))/31859534503965572279823959492121 - (921805937454099571820119166106277959316698824704*exp(-6250000000850000000001/10000000000000000000000))/121534479156362809294982755630954742431640625 + (pi^(1/2)*exp(279/16)*erfi(425000000001i/100000000000)*1053932631846001350838118399344640000i)/31859534503965572279823959492121)))");
syms y1 z
fu2=int(simplify(fu1),y1,0,z);
subplot(311);fplot(fu2,[0 10]);
subplot(312);fplot(vpa(fu2),[0 10]);
subplot(313);fplot(simplify(fu2),[0 10]);

It seems like more is happening in fplot than just a scaling issue.

[Walter Roberson]

I was thinking that perhaps when the constant was evaluated at too low of a precision, it went imaginary, leading to no plot. But no...

fu1 = str2sym("piecewise(y1 < 0, 0, 0 <= y1, -(2535301200456458802993406410752*exp(-(2*y1)/5)*exp(-(3*y1)/10)*exp(-(y1/10 + 1)^2/2)*exp(-(y1/5 + 1)^2/8)*abs(y1)^3)/(31859534503965572279823959492121*((1053932631846001350838118399344640000*pi^(1/2)*exp(279/16))/31859534503965572279823959492121 - (921805937454099571820119166106277959316698824704*exp(-6250000000850000000001/10000000000000000000000))/121534479156362809294982755630954742431640625 + (pi^(1/2)*exp(279/16)*erfi(425000000001i/100000000000)*1053932631846001350838118399344640000i)/31859534503965572279823959492121)))");
subplot(311);fplot(fu1,[0 10]);
subplot(312);fplot(real(fu1),[0 10]);
subplot(313);fplot(imag(fu1),[0 10]);

Increasing digits() to 100 does not fix the problem either.

[Paul]

Now that I look closely ... in these first two subplots it looks like there is a vertical dashed line at the edge at y1 = 10, which becomes evident with set(gca,'XLim',[0 15]). But the very first plot in this answer that you posted was the same thing, i.e., fplot(fu1 ...), but that plot went to 15, and in that plot the singularity is at y1 = 15, not y1 = 10, which also becomes evident with set(gca,...). So the location of the phantom singularity depends on the interval of fplot. There really seems to be more happneing here than just evaluation of constant terms.

>> vpa(simplify(fu1),10)
 
ans =
 
piecewise(y1 < 0, 0.0, 0 <= y1, 0.2072943186*y1^3*exp(- 0.01*y1^2 - 0.85*y1 - 0.625))

Why would fplot think there's a singularity in fu1?

[Walter Roberson]

I have had fplot show the dashed line for functions that had no singularity... or at least not when executed properly.

When I tested before with fplot(fu1) it showed a supposed singularity even if I ploted over [0 1] .

[Walter Roberson]

For this function I narrowed it down. Below 2.20150 I never saw the dashed line in the work I did. By 2.2016 I always saw it. And here with these intermediate values notice that the second has it but the third does not, indicating that it becomes dependent on exactly where you evaluate in that range.

Also, it turns out that fplot(imag(fu1)) does not show the discontinuity dash, and is always 0, so the absense of line for fplot(fu1) is not because the system goes imaginary.

fu1 = str2sym("piecewise(y1 < 0, 0, 0 <= y1, -(2535301200456458802993406410752*exp(-(2*y1)/5)*exp(-(3*y1)/10)*exp(-(y1/10 + 1)^2/2)*exp(-(y1/5 + 1)^2/8)*abs(y1)^3)/(31859534503965572279823959492121*((1053932631846001350838118399344640000*pi^(1/2)*exp(279/16))/31859534503965572279823959492121 - (921805937454099571820119166106277959316698824704*exp(-6250000000850000000001/10000000000000000000000))/121534479156362809294982755630954742431640625 + (pi^(1/2)*exp(279/16)*erfi(425000000001i/100000000000)*1053932631846001350838118399344640000i)/31859534503965572279823959492121)))");
subplot(311);fplot(fu1,[1 2.20150]);
subplot(312);fplot(fu1,[1 2.20155]);
subplot(313);fplot(fu1,[1 2.20157]);

figure
subplot(311);fplot(1./fu1,[1 2.20150]);
subplot(312);fplot(1./fu1,[1 2.20155]);
subplot(313);fplot(1./fu1,[1 2.20157]);

Ah, look at that: plotting 1/fu1 gets you a solid 0, so the implication is that fu1 is being computed as infinite!

[Christopher Creutzig]

@Paul please check https://www.mathworks.com/support/bugreports/details/2422167 .

As to the fu2 example, that one at least makes it easy to see the graph is full of numerical artefacts.

In the fu1/G case, this low-precision evaluation results in a scaling problem, because it manifests inside the scaling factor of the input:

fu1 = str2sym("piecewise(y1 < 0, 0, 0 <= y1, -(2535301200456458802993406410752*exp(-(2*y1)/5)*exp(-(3*y1)/10)*exp(-(y1/10 + 1)^2/2)*exp(-(y1/5 + 1)^2/8)*abs(y1)^3)/(31859534503965572279823959492121*((1053932631846001350838118399344640000*pi^(1/2)*exp(279/16))/31859534503965572279823959492121 - (921805937454099571820119166106277959316698824704*exp(-6250000000850000000001/10000000000000000000000))/121534479156362809294982755630954742431640625 + (pi^(1/2)*exp(279/16)*erfi(425000000001i/100000000000)*1053932631846001350838118399344640000i)/31859534503965572279823959492121)))");
G = simplify(fu1,"Steps",50)
vpa(G,5)

vpa(G)

JFTR, computing that constant in doubles gives 0.2073, and in single precision, we get 1.96e-5. The vpa(G) result is consistent with higher precision results and probably reliable.

t = @single;
t(2000000000000000000)/(t(190624999999575000000001)*exp((-t(6250000000850000000001/10000000000000000000000))) - sqrt(t(pi))*(t(831406250000000000000000)*exp(t(279/16)) - t(831406250000000000000000)*exp(t(279/16))*erf(t(425000000001/100000000000))))

[Paul]

@Christopher Creutzig

Thanks for posting the link to the bug report. However, when I run the example in the report:

>> ex = str2sym("piecewise(x < 100, (pi^2-107607675/10902937)*x)");  
fplot(ex)

I get a plot that's a horizontal line at 0, which I'm not so sure is a good example of a result that "looks reasonable." In fact, it immediately made me suspicous and I had to investigate why that was so. And then the workaround code resulted in:

ex = str2sym("piecewise(x < 100, (pi^2-107607675/10902937)*x)"); 
fplot(@(xval) subs(ex,x,xval))
Warning: Function behaves unexpectedly on array inputs. To improve performance, properly vectorize your function to return an output with the same size and shape as the input arguments. 
> In matlab.graphics.function.FunctionLine>getFunction
  In matlab.graphics.function.FunctionLine/updateFunction
  In matlab.graphics.function.FunctionLine/set.Function_I
  In matlab.graphics.function.FunctionLine/set.Function
  In matlab.graphics.function.FunctionLine
  In fplot>singleFplot (line 241)
  In fplot>@(f)singleFplot(cax,{f},limits,extraOpts,args) (line 196)
  In fplot>vectorizeFplot (line 196)
  In fplot (line 166) 
Warning: Error updating FunctionLine.

 The following error was reported evaluating the function in FunctionLine update: Undefined function or variable 'x'.

Yes, the fu2 example makes it easy to see that there is a potential problem. But the issue is whether or not the problem in the fu2 example and in the examples that @Walter Roberson posted is due to the same issue with fplot computing the constant in low precision OR if there are one or more other issues with fplot that cause spurious results (functions evaluated to inf, non-existent singularities, whatever one might call the fu2 problems) as illustrated in this thread?

[Christopher Creutzig]

@Paul The workaround section was missing a "syms x" line, thanks for checking.