# Numerical problems when plotting with fplot

1 vue (au cours des 30 derniers jours)
FRANCESCO FABRIS le 30 Mai 2023
Commenté : Paul le 30 Mai 2023
The following are two different methods to plot an exponentially modified Gaussian distribution, https://en.wikipedia.org/wiki/Exponentially_modified_Gaussian_distribution
with parameters mu, sigma, lambda.
The first one is a point-by-point calculation of the integrals based on a for cycle; the second one is direct.
WIth parameters mu=0, sigma=1, lambda=1 all works fine, and the two plots are identical. Instead, if one sets mu = 0; sigma = 2; lambda = 3.5 we have an oscillation in the second plot. If we set mu = 0; sigma = 3; lambda = 3, the second plot is completely crazy.
I cannot understand why, since we have not pointS of singularity, that can affect the numerical convergence. Note, however, that the integral between -Inf and +Inf is in any case correct, i.e. =1, since this is a probability density function.
Here are the codes:
% FIRST METHOD
syms t
mu = 0;
sigma = 2;
lambda = 3.5;
x = -5:0.1:5;
F_lambda = zeros(1,length(x));
exp_t = @(t) exp(-t.^2);
for i=1:numel(x)
xmax = +Inf;
xmin = ((mu+lambda.*sigma.^2-x(i))./(sqrt(2).*sigma));
erfc = (2./sqrt(pi)).*int(exp_t,t,xmin,xmax);
f_lambda = (lambda./2).*exp((lambda./2).*(2.*mu+lambda.*sigma.^2-2.*x(i))).*erfc;
F_lambda(i) = f_lambda;
end
plot(x,F_lambda) % SECOND METHOD
syms x t
mu = 0;
sigma = 2;
lambda = 3.5;
exp_t = @(t) exp(-t.^2);
xmax = +Inf;
xmin = ((mu+lambda.*sigma.^2-x)./(sqrt(2).*sigma));
erfc = (2./sqrt(pi)).*int(exp_t,t,xmin,xmax);
f_lambda = (lambda./2).*exp((lambda./2).*(2.*mu+lambda.*sigma.^2-2.*x)).*erfc;
double(int(@(x) f_lambda,x,-Inf,+Inf))
fplot(f_lambda) ##### 0 commentairesAfficher -1 commentaires plus anciensMasquer -1 commentaires plus anciens

Connectez-vous pour commenter.

### Réponse acceptée

Paul le 30 Mai 2023
Hi FRANCESCO,
I can't explain the inner workings of fplot, but I too have had occasion where it comes up with odd results.
Here's the second method, defining all constants as sym objects.
syms x t
mu = sym(0);
sigma = sym(2);
lambda = sym(3.5);
Define exp_t as a sym expression, didn't see a need to make it an anonymous function
exp_t = exp(-t.^2);
xmax = +Inf;
xmin = ((mu+lambda.*sigma.^2-x)./(sqrt(2).*sigma));
Define erfc expression. Matlab is smart enough to identify the integral and express erfc in term of erf
erfc = (2./sqrt(sym(pi))).*int(exp_t,t,xmin,xmax)
erfc = Here, I'll define f_lambda as a symfun.
f_lambda(x) = (lambda./2).*exp((lambda./2).*(2.*mu+lambda.*sigma.^2-2.*x)).*erfc
f_lambda(x) = %double(int(@(x) f_lambda,x,-Inf,+Inf))
Sometimes, fplot works better by increasing the MeshDensity. Didn't help here.
figure
fplot(f_lambda,'MeshDensity',200) But the Symbolic Math Toolbox has its own erfc function (and there is a numeric version as well erfc), so maybe things will work better if we use that function.
clear erfc
f_lambda(x) = (lambda./2).*exp((lambda./2).*(2.*mu+lambda.*sigma.^2-2.*x))*erfc(xmin)
f_lambda(x) = Same problem with fplot
figure
fplot(f_lambda,'MeshDensity',200) But we can plot f_lambda at specified values of x
xval = -5:0.1:5;
figure
plot(xval,double(f_lambda(xval))) ##### 2 commentairesAfficher 1 commentaire plus ancienMasquer 1 commentaire plus ancien
Paul le 30 Mai 2023
I noticed that too and it threw me for a loop. But when I used the SMT version of erfc, it came up with an expression containing -(efrc() - 2), which is the same as
(2 - erfc()) = 1 + 1 - erfc() = 1 + erf()
so I didn't pursue any further. Probably has to do with erf being an odd function so that:
erfc(a-x) = 1 - erf(a-x) = 1 + erf(x-a)

Connectez-vous pour commenter.

### Catégories

En savoir plus sur Applications dans Help Center et File Exchange

R2023a

### Community Treasure Hunt

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

Start Hunting!