Why Would fplot(f) and fplot(vpa(f)) Show Different Results?
3 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Paul
le 20 Jan 2021
Commenté : Christopher Creutzig
le 3 Fév 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 commentaires
Walter Roberson
le 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.
Réponse acceptée
Walter Roberson
le 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 commentaires
Christopher Creutzig
le 3 Fév 2021
@Paul The workaround section was missing a "syms x" line, thanks for checking.
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Assumptions dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!