MATLAB: Would fplot(f) and fplot(vpa(f)) Show Different Results

fplotMATLABSymbolic Math Toolbox

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))

Best Answer

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.
Related Question