Numerical Integration Is my approach correct?

I have an integration which can not be solved in closed form and I have to solve numerically. I am using MATLAB built functions integrate. My function which I am trying to integrate over is in function of 3 parameters ( v, gamma1,gamma2) it is is a one integral with dummy variable v, my integration limit is from gamma2 to infinity.

Should I run a for loop over gamma1 and gamma2 ? That is should I do the following, ( note that v is my dummy variable, and my function has gamma1 gamma2 and v$

    gamma1=1:1:1000;
    gamma2=1:1:1000;
    for i=1:1:length(gamma1)
    for j=1:1:length(gamma2)
    gamma12(i,j)= gamma1(i)^(-1)+gamma2(j)^(-1);  % part of my function
    g(i,j)= v^(-1)*(10/gamma12(i,j));            % part of my function
    fun(i,j)= exp(v)*(g(i,j)/(1+g(i,j))) ;% the function i want to integ
    y1{i,j}=matlabFunction(fun(i,j));
    y(i,j)=integrate (fun(i,j),gamma2(j),+inf)
    end
    end
I guess my question is, do you think what I do makes sense?

thanks

Answer 1

[Mike Hosea]

Well, on first reading, I just wanted to make what you wrote work, and you will find some help below on that, but on the second reading, I wonder why you are looping over so many values of gamma1 and gamma2. What are you going to do with all this data? You will be creating a 1000-by-1000 matrix. This could take awhile. I suggest that you make it 10-by-10 temporarily until you are sure that it is working! Then you can define gamma1 and gamma2 to be whatever you happen to need.

Anyway, here is what you need to make the loops work.

Instead of defining the large matrix gamma12(i,j) in the loop, you probably want to just create a gamma12 scalar that changes in each iteration. You could write

for i=1:1:length(gamma1)
    for j=1:1:length(gamma2)
        gamma12 = gamma1(i)^(-1) + gamma2(j)^(-1);
        <snip>

Now that brings us to the definition of a function. You should not need the matlabFunction function. Instead, define your functions using the @ operator. When you use @ followed by a list of variables (comma separated, with bracketing parentheses), those variables are considered arguments. For example, @(x,y,z) means that we are defining a function of variables x, y, and z. Here we just have one variable, v, so we will use @(v) to indicate that what follows is to be regarded as a function of v. So your g function would look like this

g = @(v)v^(-1)*(10/gamma12);

But I can tell you right now that we're going to need to make a small change. The reason is that v^(-1) is matrix inversion. But what you meant was reciprocation. They're the same thing for scalars, very different for matrices, and the INTEGRAL function doesn't just evaluate one element at a time--it passes arrays of elements that it expects you to evaluate in a batch, element-by-element. So we just need to change this to .^, the element-wise power function.

g = @(v)v.^(-1)*(10/gamma12);

Now we can define

fun = @(v)exp(v).*(g(v)./(1 + g(v)));

Notice how I used .* multiplication instead of * multiplication. That's because * is matrix multiplication, not element-by-element multiplication, and I used ./ instead of / because / is matrix "division" (solving linear systems), whereas ./ is element-by-element division.

Now, you'll be ready to define

y(i,j) = integral(fun,gamma2,+inf)

Note that it is "integral", not "integrate". See if you can make that work.

[DM]

Thanks for your detailed reply, after reading your suggestions, I tried the following code

gamma1=1:1:10;
gamma2=1:1:10;

    for i=1:length(gamma1)

      for j=1:length(gamma2)

gamma12 = gamma1(i)^(-1) + gamma2(j)^(-1);
g = @(v)v.^(-1)*(10/gamma12);
fun = @(v)exp(v).*(g(v)./(1 + g(v)));
y(i,j) = integral(fun,gamma2,+inf)

             end
      end

 y(i,j) = integral(fun,gamma2(j),+inf)

Do you think I did the right thing?

[Mike Hosea]

Well, you have an extra y(i,j) = integral... at the bottom to delete, and the y(i,j) = integral line inside the loop should read

y(i,j) = integral(fun,gamma2(j),+inf);

(you need to subscript gamma2--sorry for messing that up before--I changed my mind about the loop variables and forgot to re-edit that).

Also, surely exp(v) should be exp(-v). Otherwise, the integral won't converge.

After that, I think it will run, but whether it is the right thing to do depends on what you are trying to accomplish. When all is said and done, you will have a 10-by-10 array of numbers, each corresponding to a particular gamma1(i) and gamma2(j) pair. Hopefully you know what you want to do with those numbers.