Help on double integral

49 views (last 30 days)
Mohamed on 31 Aug 2014
Hello Everyone,
I need your help on this problem. I am new to Matlab and still building my skills. Here is issue. I want to double integrate the following:
theta limits = [0,pi/2] phi limits = [0,2*pi]
k is a variable that ranges from N:M
function = cos(k.*sin(theta).*cos(phi))
I want to double integrate the above function with respect to theta and phi limits. I need to repeat the integration every time k change its value based on the k range above.
dblquad works if k has one value but I couldn't use it with k range. Any help would be highly appreciated.
  1 Comment
Mohamed on 31 Aug 2014
Just to add to my question. I am using an old Matlab software Ver 6 Release 12

Sign in to comment.

Answers (2)

Star Strider
Star Strider on 31 Aug 2014
This works:
fun = @(k,theta,phi) cos(k.*sin(theta).*cos(phi));
kv = 5:8;
for k1 = 1:length(kv)
k = kv(k1);
ifun(k1,:) = [k dblquad(@(theta,phi) fun(k,theta,phi), 0, 2*pi, 0, 2*pi)];
This defines k initially as kv (here, arbitrarily from 5 to 8), then loops through every value of k. The ‘ifun’ matrix stores the value of k and the value of the integral in each row.
The quad2d function is faster than dblquad, and appears to give the same results.
Star Strider
Star Strider on 31 Aug 2014
First, if you don’t want to store k in the result, the ‘ifun’ assignment line becomes:
ifun(k1) = dblquad(@(theta,phi) fun(k,theta,phi), 0, 2*pi, 0, 2*pi);
I have no idea what that error is. I’m using anonymous functions, and that code worked perfectly for me. (If it hadn’t, I’d not have posted it.) The complete function should be entered as:
fun = @(k,theta,phi) cos(k.*sin(theta).*cos(phi));
I don’t remember when anonymous functions were introduced. If the anonymous function doesn’t work in your version, you will have to create a function file:
function val = fun(k,theta,phi)
val = cos(k.*sin(theta).*cos(phi));
then save it as fun.m (or rename the function to something less ambiguous and save it to the same file name as the function name).
You also may have to remove k as a function argument and substitute it manually into the function (or take the risk of declaring k as a global variable inside and outside the function for the purposes of this integration). The function file would then be:
function val = fun(theta,phi)
global k
val = cos(k.*sin(theta).*cos(phi));
I tested the integration with an inline function and couldn’t get it to work. If you cannot use an anonymous function, you will have to use a function file to do the integration, and a global variable for k. I don’t have access to R12, so I can’t test my code with it. Unfortunately, you will have to experiment. With luck, one of my alternative suggestions will work for you.

Sign in to comment.

Mike Hosea
Mike Hosea on 26 Sep 2014
Edited: Mike Hosea on 26 Sep 2014
ffull = @(k,theta,phi)cos(k.*sin(theta).*cos(phi));
qscalar = @(kscalar)integral2(@(theta,phi)ffull(kscalar,theta,phi),0,pi/2,0,2*pi);
q = @(k)arrayfun(qscalar,k);
Now what do you want to do with the q function?
>> q(1:3)
ans =
8.692413142202774 5.778924787670234 2.585516200656736
or maybe plot it for k = 1:20
>> k = 1:20;
>> plot(k,q(k))
A key technical point that I've sort of glossed over here is that it is critical that you "bind" the k value by creating a function handle after a particular, scalar value of k has been selected. If you loop over the k values you need, you would create the function of just two variables, theta and phi, inside the loop. Here it was convenient to create it inside the call to INTEGRAL2. Another technical point is that because INTEGRAL2 can only calculate one value at a time, you have to use ARRAYFUN to vectorize it. You don't have to vectorize it if instead you decide to call it inside a loop over the k values.
  1 Comment
José Anchieta de Jesus Filho
so, I need to create a loop that calculates integrals at each point. Say I have the following function handle
D = @(x,y) x.*(Z(x,y).^2);
from that function, I need to calculate a double integral and create the following loop
int1 = @(y) integral(@(x)D(x,y),x1,x2,'ArrayValued',true);
int2 = integral(@(y)int1,y1,y2,'ArrayValued',true);
c1 = int2./g % g is a scalar number
int11 = @(y) integral(@(x)D(x,y),x2,x3,'ArrayValued',true);
int22 = integral(@(y)int11,y1,y2,'ArrayValued',true);
c11 = int22./g % g is a scalar number
int111 = @(y) integral(@(x)D(x,y),x3,x4,'ArrayValued',true);
int222 = integral(@(y)int111,y1,y2,'ArrayValued',true);
c111 = int222./g % g is a scalar number
i'm trying to do like this
x = 1:0.1:5
for i = 1:length(x)
D = @(x,y) x.*(Z(x,y).^2);
int1(i) = @(y) integral(@(x)D(x,y),x(i),x(i)+0.5,'ArrayValued',true);
int2(i) = integral(@(y)int1(i),y1,y2,'ArrayValued',true);
c = int2(i)./g
I need to plot a graph of c versus x

Sign in to comment.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!