Numerical Integration Is my approach correct?

1 view (last 30 days)
DM
DM on 27 Oct 2014
Edited: Mike Hosea on 29 Oct 2014
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

Answers (1)

Mike Hosea
Mike Hosea on 28 Oct 2014
Edited: Mike Hosea on 28 Oct 2014
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.
  2 Comments
Mike Hosea
Mike Hosea on 29 Oct 2014
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.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!