Why Would fplot(f) and fplot(vpa(f)) Show Different Results?
4 views (last 30 days)
Show older comments
Paul
on 20 Jan 2021
Commented: Christopher Creutzig
on 3 Feb 2021
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
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.
Accepted Answer
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)))")
G = simplify(fu1,'steps', 50)
%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))
vpa(subs(fu1,y1,1))
double(subs(vpa(fu1),y1,1))
%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))
vpa(subs(G,y1,1))
double(subs(vpa(G),y1,1))
%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
Christopher Creutzig
on 3 Feb 2021
@Paul The workaround section was missing a "syms x" line, thanks for checking.
More Answers (0)
See Also
Categories
Find more on Assumptions in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!