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

4 views (last 30 days)
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))
  2 Comments
Walter Roberson
Walter Roberson on 20 Jan 2021
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
Paul on 20 Jan 2021
Edited: Paul on 20 Jan 2021
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

Sign in to comment.

Accepted Answer

Walter Roberson
Walter Roberson on 20 Jan 2021
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.
  20 Comments
Paul
Paul on 3 Feb 2021
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
Christopher Creutzig on 3 Feb 2021
@Paul The workaround section was missing a "syms x" line, thanks for checking.

Sign in to comment.

More Answers (0)

Tags

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!