Evaluate integral of vectorial function

Dear all i want to evaluate an integral for a function with different values for some parameters but, instead of performing it in a loop, i would like to do it vectorializing

fun = @(x,c) 1./(x.^3-2*x-c);
q = integral(@(x)fun(x,c),0,c)

where the parameter c is a mx1 vector. and the result q of the integration should be a mx1 vector as well. However, if i define c as a vector the integration does not work. In mycase, performing a the integral into a loop is very time consuming. Is there another way?

[John BG]

are you going to accept an answer among the ones supplied?

would you please be so kind to comment?

Answer 1

[John BG]

if you choose c values that include or are too close to a pole, the function integral may trigger error warning because either from the method used or because the integral does not exist the error limit has reached 0.21

what a bout the following:

fun = @(x,c) 1./(x.^3-2*x-c);
c=[1:10]';

qv=[];
for k=1:1:numel(c) 
   q=integral(@(x)fun(x,c(k)),0,c(k));
   qv=[qv q]; 
end

qv
=
  Columns 1 through 4
  -0.598079642332067  -0.122891186512513   5.687368137358827  -0.304644136939809
  Columns 5 through 8
  -0.263556686781940  -1.859012022895603   0.688146342208949   0.052739282210456
  Columns 9 through 10
  -0.069075724263797  -0.236803884145312

each loop triggers the following warning

Warning: Reached the limit on the maximum number of intervals in use. Approximate
bound on error is   2.1e-01. The integral may not exist, or it may be difficult to
approximate numerically to the requested accuracy. 
> In integralCalc/iterateScalarValued (line 372)
  In integralCalc/vadapt (line 132)
  In integralCalc (line 75)
  In integral (line 88) 

 for c=-1

roots([1 0 -2 -1])
ans =
     1.618033988749895
    -1.000000000000001
    -0.618033988749895

if any of the c values is too close one of the poles the function integral may return the warning kind of telling not really sure the returned value is really the integral or if even exists

the location of the poles as function of c

syms c
>> roots([1 0 -2 c])
ans =
                                                                                                      2/(3*((c^2/4 - 8/27)^(1/2) - c/2)^(1/3)) + ((c^2/4 - 8/27)^(1/2) - c/2)^(1/3)
 - 1/(3*((c^2/4 - 8/27)^(1/2) - c/2)^(1/3)) - ((c^2/4 - 8/27)^(1/2) - c/2)^(1/3)/2 - (3^(1/2)*(2/(3*((c^2/4 - 8/27)^(1/2) - c/2)^(1/3)) - ((c^2/4 - 8/27)^(1/2) - c/2)^(1/3))*1i)/2
 - 1/(3*((c^2/4 - 8/27)^(1/2) - c/2)^(1/3)) - ((c^2/4 - 8/27)^(1/2) - c/2)^(1/3)/2 + (3^(1/2)*(2/(3*((c^2/4 - 8/27)^(1/2) - c/2)^(1/3)) - ((c^2/4 - 8/27)^(1/2) - c/2)^(1/3))*1i)/2

the function integral [-Inf,Inf] varies with the position of the poles

qv=[];
for k=-10:.01:10
fun = @(x) 1./(x.^3-2*x-k);
q=integral(fun,-Inf,Inf);
qv=[qv q];
end
plot(qv)

.

I have attached to this answer the variable qv in file qv.mat because it takes a while to generate [-Inf,Inf] integral values with c as parameter between [-10 10] with c resolution 0.01

So may be you would like to consider integrating fun between [0 c] by doing

qv(c1)-qv(0)

if you find my answer useful would you please mark it as Accepted Answer by clicking on the ACCEPT ANSWER button?

thanks in advance for time and attention

John BG

<mailto:[e-mail] [e-mail]>

Answer 2

[Star Strider]

The arrayfun function may do what you want.

See if this works for you:

c = 1:5;                                                    % Create ‘c’
fun = @(x,c) 1./(x.^3-2*x-c);
q = @(c) integral(@(x)fun(x,c),0,c);
r = arrayfun(q,c)

Answer 3

[David Goodmanson]

Hello Davide, I believe an explicit solution is

cvec = .01:.01:20;
q = []; 
for c = cvec
R = roots([1 0 -2 -c]);
r = R(1);  s = R(2);  t = R(3);
I = real( ((r-s)*(r-t))^-1*log((c-r)/(-r)) ...
         +((s-t)*(s-r))^-1*log((c-s)/(-s)) ...
         +((t-r)*(t-s))^-1*log((c-t)/(-t)) );
q = [q I];         % should really preallocate instead
end
plot(cvec,q);

There is a pole in q at c = sqrt(3). The values of q to the right of sqrt(3) are probably correct if q is taken to be the principal value of the integral. The code does not work for c exactly at zero.

Answer 4

[Walter Roberson]

Your question is ambiguous. Are you doing the equivalent of

   for K = 1 : length(c)
     this_c = c(K);
     result(K) = integral( @(x) fun(x, this_c), 0, this_c);
   end

where you are producing one output for each c value, with the individual value of c used in the calculation and also as the upper bound for the same calculation?

Or are you doing the equivalent of

   for K_outer = 1 : length(c)
     c_outer = c(K_outer);
     for K_inner = 1 : length(c)
       c_inner = c(K_inner)
       result(K_outer, K_inner) = integral( @(x) fun(x, c_inner), 0, c_outer);
     end
   end

where you are producing one output for each combination of c values?

MATLAB's integral() function cannot handle either possibility with a single step. It can, however, do

   for K_outer = 1 : length(c)
     c_outer = c(K_outer);
     result(K_outer, :) = integral( @(x) fun(x, c), 0, c_outer, 'ArrayValued', true);
   end

This unfortunately cannot be converted to handle the case where each corresponding value is used in the function and as the upper bound.