Evaluate integral of vectorial function
2 visualizaciones (últimos 30 días)
Mostrar comentarios más antiguos
Davide Leonetti
el 25 de Dic. de 2016
Comentada: John BG
el 28 de Dic. de 2016
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?
1 comentario
John BG
el 28 de Dic. de 2016
are you going to accept an answer among the ones supplied?
would you please be so kind to comment?
Respuesta aceptada
John BG
el 25 de Dic. de 2016
Editada: John BG
el 25 de Dic. de 2016
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)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/176258/image.jpeg)
.
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
0 comentarios
Más respuestas (3)
Star Strider
el 25 de Dic. de 2016
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)
0 comentarios
David Goodmanson
el 26 de Dic. de 2016
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.
0 comentarios
Walter Roberson
el 27 de Dic. de 2016
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.
0 comentarios
Ver también
Categorías
Más información sobre Symbolic Math Toolbox en Help Center y File Exchange.
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!