Thread Subject:

Problem with results from symbolic toolbox

Here is the entirety of the offending code:

function legendre_plot(lambda,delta)

xsize=9;
xtest= (1 ./ (xsize+1)) .* (1:xsize);
legendre_data = zeros(1,xsize);

for j=1:xsize
    legendre_data(j) = value_legendre_1(xtest(j));
end

plot(xtest,legendre_data)

function out = value_legendre_1(x)
t = sym('t' ,'real');
legendre_t = solve(diff(value_cgf_1(t),t)-x);
out = real(legendre_t .* x - value_cgf_1(legendre_t));
end

function out = value_cgf_1(t)
y = sym('y' ,'real');
out =log(int(exp(t.*y)*value_pdf_1(y),y,0,1));
end

function out = value_pdf_1(x)
out = -((lambda .* x.^(-((lambda + log(delta))./log(delta))))./
log(delta));
end

end


When I run this for lambda=(2/7)*(1/200), delta = .95, two values I am
interested in, I get a graph that looks about how I think it should.
When I run this for for lambda = (2/7)(1-(1/200)), delta = .95, two
other values I am interested in, I get the following error:

??? The following error occurred converting from sym to double:
Error using ==> mupadmex
Error in MuPAD command: DOUBLE cannot convert the input expression
into a
double array.

If the input expression contains a symbolic variable, use the VPA
function
instead.

Error in ==> legendre_plot at 17
    legendre_data(j) = value_legendre_1(xtest(j));


What is going on here?

On 07.09.11 18:37, DOD wrote:
> What is going on here?

The int call did not return something in closed form; the result still
contains a symbolic “limit” expression. In this case, you are lucky; the
problematic integration is one you can replace with numerical quadrature:

function legendre_plot(lambda,delta)

xsize=9;
xtest= (1 ./ (xsize+1)) .* (1:xsize);
legendre_data = zeros(1,xsize);

for j=1:xsize
    legendre_data(j) = value_legendre_1(xtest(j));
end

plot(xtest,legendre_data)

function out = value_legendre_1(x)
t = sym('t' ,'real');
legendre_t = double(solve(diff(value_cgf_1(t),t)-x));
out = real(legendre_t .* x - ...
  quad(@(y) exp(legendre_t.*y).*double(value_pdf_1(y)),0,1));
end

function out = value_cgf_1(t)
y = sym('y' ,'real');
out =log(int(exp(t.*y)*value_pdf_1(y),y,0,1));
end

function out = value_pdf_1(x)
out = -((lambda .* x.^(-((lambda + log(delta))./ ...
         log(delta))))./ log(delta));
end

end



Christopher

On Sep 8, 7:46 am, Christopher Creutzig
<Christopher.Creut...@mathworks.com> wrote:
> On 07.09.11 18:37, DOD wrote:
>
> > What is going on here?
>
> The int call did not return something in closed form; the result still
> contains a symbolic “limit” expression. In this case, you are lucky; the
> problematic integration is one you can replace with numerical quadrature:
>
> function legendre_plot(lambda,delta)
>
> xsize=9;
> xtest= (1 ./ (xsize+1)) .* (1:xsize);
> legendre_data = zeros(1,xsize);
>
> for j=1:xsize
>     legendre_data(j) = value_legendre_1(xtest(j));
> end
>
> plot(xtest,legendre_data)
>
> function out = value_legendre_1(x)
> t = sym('t' ,'real');
> legendre_t = double(solve(diff(value_cgf_1(t),t)-x));
> out = real(legendre_t .* x - ...
>   quad(@(y) exp(legendre_t.*y).*double(value_pdf_1(y)),0,1));
> end
>
> function out = value_cgf_1(t)
> y = sym('y' ,'real');
> out =log(int(exp(t.*y)*value_pdf_1(y),y,0,1));
> end
>
> function out = value_pdf_1(x)
> out = -((lambda .* x.^(-((lambda + log(delta))./ ...
>          log(delta))))./ log(delta));
> end
>
> end
>
> Christopher

You definitely found the problem, the integral is not as nice as I had
hoped. the numerical approach doesn't quite solve my problem, as this
new program illustrates:

function [symbolic_version quadl_version quadgk_version quad_version]
= problematic_integrals
z = sym('z', 'real');
t = -0.8566;
lambda = 1/700;
delta = .95;

symbolic_version = double(int(exp(t.*z).*value_pdf_1(z,lambda, delta),
0,1));
quadl_version = quadl(@(x)exp(t.*x) .* value_pdf_1(x,lambda,delta),
0,1);
quadgk_version = quadgk(@(x)exp(t.*x) .* value_pdf_1(x,lambda,delta),
0,1);
quad_version = quad(@(x)exp(t.*x) .* value_pdf_1(x,lambda,delta),0,1);


    function out = value_pdf_1(x,lambda, delta)

    out = -((lambda .* x.^(-((lambda + log(delta))./log(delta))))./
log(delta));
    end

end

I get a different answer from each program. The one I believe to be
correct is int. What can I do to get correct answers out of a
numerical version? Even with a singularity at 0, I figured this
integral is basically pretty well behaved, but I guess not.

Thanks,

Dennis

On Sep 8, 4:44 pm, DOD <dco...@gmail.com> wrote:
> On Sep 8, 7:46 am, Christopher Creutzig
>
>
>
>
>
>
>
>
>
> <Christopher.Creut...@mathworks.com> wrote:
> > On 07.09.11 18:37, DOD wrote:
>
> > > What is going on here?
>
> > The int call did not return something in closed form; the result still
> > contains a symbolic “limit” expression. In this case, you are lucky; the
> > problematic integration is one you can replace with numerical quadrature:
>
> > function legendre_plot(lambda,delta)
>
> > xsize=9;
> > xtest= (1 ./ (xsize+1)) .* (1:xsize);
> > legendre_data = zeros(1,xsize);
>
> > for j=1:xsize
> >     legendre_data(j) = value_legendre_1(xtest(j));
> > end
>
> > plot(xtest,legendre_data)
>
> > function out = value_legendre_1(x)
> > t = sym('t' ,'real');
> > legendre_t = double(solve(diff(value_cgf_1(t),t)-x));
> > out = real(legendre_t .* x - ...
> >   quad(@(y) exp(legendre_t.*y).*double(value_pdf_1(y)),0,1));
> > end
>
> > function out = value_cgf_1(t)
> > y = sym('y' ,'real');
> > out =log(int(exp(t.*y)*value_pdf_1(y),y,0,1));
> > end
>
> > function out = value_pdf_1(x)
> > out = -((lambda .* x.^(-((lambda + log(delta))./ ...
> >          log(delta))))./ log(delta));
> > end
>
> > end
>
> > Christopher
>
> You definitely found the problem, the integral is not as nice as I had
> hoped.  the numerical approach doesn't quite solve my problem, as this
> new program illustrates:
>
> function [symbolic_version quadl_version quadgk_version quad_version]
> = problematic_integrals
> z = sym('z', 'real');
> t = -0.8566;
> lambda = 1/700;
> delta = .95;
>
> symbolic_version = double(int(exp(t.*z).*value_pdf_1(z,lambda, delta),
> 0,1));
> quadl_version = quadl(@(x)exp(t.*x) .* value_pdf_1(x,lambda,delta),
> 0,1);
> quadgk_version = quadgk(@(x)exp(t.*x) .* value_pdf_1(x,lambda,delta),
> 0,1);
> quad_version = quad(@(x)exp(t.*x) .* value_pdf_1(x,lambda,delta),0,1);
>
>     function out = value_pdf_1(x,lambda, delta)
>
>     out = -((lambda .* x.^(-((lambda + log(delta))./log(delta))))./
> log(delta));
>     end
>
> end
>
> I get a different answer from each program.  The one I believe to be
> correct is int.  What can I do to get correct answers out of a
> numerical version?  Even with a singularity at 0, I figured this
> integral is basically pretty well behaved, but I guess not.
>
> Thanks,
>
> Dennis

I have realized that these integrals can be presented by a combination
of gamma functions. For anyone's future reference, this code gives
the exact values for the log of those integrals(which is what I want
anyway):


function out = new_exact_cgf(t,lambda,delta)

A = - lambda ./ log(delta);

out = real(log( A .* ((-t) .^ (-A)) .* (gamma(A) - gamma(A).*gammainc(-
t,A,'upper'))));

end

I suppose the whole reason these integral are implemented as special
functions is the numerical routines have trouble with them.

Dennis