[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: Numerical integration with quadcc (resubmitted with MWE)

From: Przemek Klosowski
Subject: Re: Numerical integration with quadcc (resubmitted with MWE)
Date: Fri, 7 Jun 2019 18:21:02 -0400
User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:60.0) Gecko/20100101 Thunderbird/60.2.1

On 6/7/19 5:14 PM, BGreen wrote:
My problem is that I get the error "integrand F must return a single,
real-valued vector of the same size as the input" even when the integrand
DOES return such a vector. I am using quadcc because I have infinite limits
of integration and a piecewise-defined integrand. Things will execute
correctly if I just use quad, but I'm not sure whether or not the answer
will be trustworthy since the documentation recommends other integrators for
infinite bounds and nonsmooth functions.

I think the problem results from the fact that your integrand returns a NxN matrix when given a length-N column vector.

I added a printout of a test case before quadcc():

function Fsame = IntDemonstration(lB,edielec,d,D)
>         Vsame = @(x) (1./edielec).*Vsame_fit(x./lB);
>         F0011 = @(x) (1-(x.^2)/2).*exp(-(x.^2./2));
>         integ_F0011same = @(x) F0011(x).*Vsame(x);
>         integ_F0011same ([1;2])
>         Fsame(1,1,2,2) = quadcc(integ_F0011same,0.000001,Inf);
> end
octave:80> IntDemonstration(256.56/sqrt(25),1,3.4,200)
ans =

   0.29696   0.29299
  -0.13252  -0.13075

I am sorry but I ran out of time to figure out where this happens---you seem to carefully do element by element operations, so I must be missing one place where the auto-broadcast happens...

I call the code below with "Fsame =
IntDemonstration(256.56/sqrt(25),1,3.4,200);". Everything below is in the
same file, "IntDemonstration.m".


function Fsame = IntDemonstration(lB,edielec,d,D)
        Vsame = @(x) (1./edielec).*Vsame_fit(x./lB);
        F0011 = @(x) (1-(x.^2)/2).*exp(-(x.^2./2));
        integ_F0011same = @(x) F0011(x).*Vsame(x);
        Fsame(1,1,2,2) = quadcc(integ_F0011same,0.000001,Inf);

function y = Vsame_fit(q)
        for k=1:length(q)
                if q(k)<0.022913
                        y(k) = 0.98*tanh(200*q(k));
                elseif q(k)<0.50274
                        y(k) = 1/(0.9*q(k)+1);
                elseif q(k)<21.598
                        y(k) = 1/(1.046*q(k)+0.9266);
                        y(k) = 1/(0.9512*q(k)+2.89);


I wrote the loop over the length of the argument to ensure that it would
return a vector. Furthermore, I tried inserting something simple as a test.
I put the following code in right after the definition of integ_F0011same:


As expected, Octave prints out the 10-element vector that results from that
computation. Then it throws an error when trying to evaluate the next line
with quadcc.

At the end of the day, I just want to have the value of the integral which I
was hoping would be saved to Fsame(1,1,2,2). For the record, the reason I
have defined things with so many separate functions is that I reuse Vsame
three more times to evaluate similar integrals, e.g.

integ_F0110same = @(x) F0110(x).*Vsame(x);
Fsame(1,2,2,1) = quadcc(integ_F0110same,0,Inf);

where F011(x) was also defined with a function handle. Is there a better way
I should be doing this?

Sent from:

reply via email to

[Prev in Thread] Current Thread [Next in Thread]