Discover MakerZone

MATLAB and Simulink resources for Arduino, LEGO, and Raspberry Pi

Learn more

Discover what MATLAB® can do for your career.

Opportunities for recent engineering grads.

Apply Today

Thread Subject:
Problem with results from symbolic toolbox

Subject: Problem with results from symbolic toolbox

From: DOD

Date: 7 Sep, 2011 16:37:38

Message: 1 of 4

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?

Subject: Problem with results from symbolic toolbox

From: Christopher Creutzig

Date: 8 Sep, 2011 12:46:48

Message: 2 of 4

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

Subject: Problem with results from symbolic toolbox

From: DOD

Date: 8 Sep, 2011 21:44:37

Message: 3 of 4

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

Subject: Problem with results from symbolic toolbox

From: DOD

Date: 9 Sep, 2011 19:45:46

Message: 4 of 4

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

Tags for this Thread

No tags are associated with this thread.

What are tags?

A tag is like a keyword or category label associated with each thread. Tags make it easier for you to find threads of interest.

Anyone can tag a thread. Tags are public and visible to everyone.

Contact us