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

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.

Christopher Creutzig
on 3 Feb 2021 at 7:55

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

Christopher Creutzig
on 3 Feb 2021 at 14:17

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

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

Start Hunting!
## 2 Comments

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/721739-why-would-fplot-f-and-fplot-vpa-f-show-different-results#comment_1273684

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/721739-why-would-fplot-f-and-fplot-vpa-f-show-different-results#comment_1273684

## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/721739-why-would-fplot-f-and-fplot-vpa-f-show-different-results#comment_1273699

⋮## Direct link to this comment

https://au.mathworks.com/matlabcentral/answers/721739-why-would-fplot-f-and-fplot-vpa-f-show-different-results#comment_1273699

Sign in to comment.