# Double integration of a binomial function (in x and y) generated from two row arrays

8 views (last 30 days)
Angshuman Sharma on 27 Jun 2020
Commented: John D'Errico on 29 Jun 2020
I have two row arrays containing two different variables, e.g. ax is a row array expressed as, ax = x - a, where x is a variable and a is a constant.
by is the other row array expressed as, by = y - b, where y is a variable and b is a constant.
I have mathematically processed these two row arrays to generate an algebraic binomial function (which is not an array) in x and y.
Now, I need to perform definite double integration of this polynomial in x and y over 0 < x < A, and 0 < y < B.
How can I do that? Any help will be much appreciated.
Note: I tried using integral2(fun, 0, A, 0, B), but it did not work, because the binomial function generated by MATLAB contains algebraic terms that are multiplied (*), divided (/), etc., e.g. P = (x-a1)*(x-a2)/(x-a3)
Now, integral2 needs me to use (.) before *,/, ^ operations, like P = (x-a1).*(x-a2)./(x-a3), which is a hurdle to do manually as the expression is huge.
I want MALAB to do the integration without my intervention.

darova on 27 Jun 2020
Please show the code you tried. Show the formula
Angshuman Sharma on 27 Jun 2020
Thank you all for trying to help me with my problem. The code that I am using is as follows:
A = 10;
B = 20;
a = [0, 1, 2, 3, 4];
b = [5, 6, 7, 8, 9];
syms x y
x = sym('x');
y = sym('y');
ax = x - a;
by = y - b;
r = ax.^2 + by.^2
K1 = - by./(2*pi.*r);
K2 = ax./(2*pi.*r);
K1sum = sum(K1, 'all');
K2sum = sum(K2, 'all');
Ksq = K1sum^2 + K2sum^2;
K = 0.5*Ksq;
syms f(x, y)
f(x, y) = @(x, y) K;
P = integral2(Ep, -A, A, -B, B)
I need to integrate K, which is a function of x and y, over their respective limits. I am so bad with coding.
When I run this code, it gives me the following error:
Error using integral2 (line 82)
First input argument must be a function handle.
Error in integration (line 53)
P = integral2(Ep, -A, A, -B, B)

John D'Errico on 27 Jun 2020
I follow what you did (possibly)
A = 10;
B = 20;
a = [0, 1, 2, 3, 4];
b = [5, 6, 7, 8, 9];
syms x y
x = sym('x');
y = sym('y');
ax = x - a;
by = y - b;
r = ax.^2 + by.^2
K1 = - by./(2*pi.*r);
K2 = ax./(2*pi.*r);
up to here. K1 and K2 are:
K1
K1 =
[ -(y - 5)/(2*pi*((y - 5)^2 + x^2)), -(y - 6)/(2*pi*((x - 1)^2 + (y - 6)^2)), -(y - 7)/(2*pi*((x - 2)^2 + (y - 7)^2)), -(y - 8)
K2
K2 =
[ x/(2*pi*((y - 5)^2 + x^2)), (x - 1)/(2*pi*((x - 1)^2 + (y - 6)^2)), (x - 2)/(2*pi*((x - 2)^2 + (y - 7)^2)), (x - 3)/(2*pi*((x - 3)^2 + (y - 8)^2)), (x - 4)/(2*pi*((x - 4)^2 + (y - 9)^2))]
So K1 and K2 are VECTORS. Next, you create K1sum.
K1sum = sum(K1, 'all');
Why do you think you need to use the all option? The reason I ask is because when you do something strange, that suggests you may not know what you are doing, or why you did it. It also suggests that you may be expecting something completely different.
And since I have been given no clue as to what you are doing here, or why you are doing it, I can only look at your code and wonder whether you do not understand sum, or you are expecting something else to have happened. Since sum, when applied to a vector, already sums all elements of that vector, there is no need to use the 'all' option. So was this just an overly zealous use of sum?
Assuming that you do know what you are doing (should I?)
K1sum = sum(K1, 'all');
K2sum = sum(K2, 'all');
Ksq = K1sum^2 + K2sum^2;
K = 0.5*Ksq;
So K is a symbolic scalar variable.
Next, you convert K into a function, so you can use integral2. Part of me wonders why you are using integral2, instead of just calling int twice to integrat ove the rectangle of interest. But then I look at K, which is to be polite, a bit of a mess. Were I to expand that in an attept to integrate it, I'd expect int to take a long time to return.
pretty(K)
/ x - 1 x - 2 x - 3 x - 4 x \2 / y - 6 y - 7 y - 8 y - 9 y - 5 \2
| ----- + ----- + ----- + ----- + -- | | ----- + ----- + ----- + ----- + ----- |
\ #4 #3 #2 #1 #5 / \ #4 #3 #2 #1 #5 /
--------------------------------------- + ------------------------------------------
2 2
where
2 2
#1 == 2 pi ((x - 4) + (y - 9) )
2 2
#2 == 2 pi ((x - 3) + (y - 8) )
2 2
#3 == 2 pi ((x - 2) + (y - 7) )
2 2
#4 == 2 pi ((x - 1) + (y - 6) )
2 2
#5 == 2 pi ((y - 5) + x )
So integral2 it is.
This means, IF I can assume you knoew what you were doing, before, is you just have no idea how to convery K into a function that integral2 can integrate.
Kxy = matlabFunction(K)
Kxy =
function_handle with value:
@(x,y)((x-1.0)./(pi.*((x-1.0).^2+(y-6.0).^2).*2.0)+(x-2.0)./(pi.*((x-2.0).^2+(y-7.0).^2).*2.0)+(x-3.0)./(pi.*((x-3.0).^2+(y-8.0).^2).*2.0)+(x-4.0)./(pi.*((x-4.0).^2+(y-9.0).^2).*2.0)+x./(pi.*((y-5.0).^2+x.^2).*2.0)).^2./2.0+((y-6.0)./(pi.*((x-1.0).^2+(y-6.0).^2).*2.0)+(y-7.0)./(pi.*((x-2.0).^2+(y-7.0).^2).*2.0)+(y-8.0)./(pi.*((x-3.0).^2+(y-8.0).^2).*2.0)+(y-9.0)./(pi.*((x-4.0).^2+(y-9.0).^2).*2.0)+(y-5.0)./(pi.*((y-5.0).^2+x.^2).*2.0)).^2./2.0
Kxy(2,3)
ans =
0.017393
As you see, I can evaluate it at any combination of x and y.
But before I throw that into integral2, PLOT IT!
fsurf(Kxy,[-A,A -B,B]) So almost entirely zero, except for a series of delta-like spikes. Sigh. Do I really want to take the next step?
integral2(Kxy,-A,A,-B,B)
Warning: Non-finite result. The integration was unsuccessful. Singularity likely.
> In integral2Calc>integral2t (line 121)
In integral2Calc (line 9)
In integral2 (line 106)
ans =
NaN
Why did integral2 fail? Because those spikes are a bit nasty. Not only is there a singularity, but it involves something like 0/0 at the center of those singularities.
Kxy(0,5)
ans =
NaN
Sorry, but I'm starting to wonder if this is something you really wanted to compute, or if you have made mistakes, etc. In that case, I might just be wasting a fair amount of time trying to resolve those singularities. The question would need to be asked if the integral exists at all.

Angshuman Sharma on 27 Jun 2020
What I showed you was just an example that I created taking some random values for a, b, A and B. This is only a tiny part of a very long code (like the tip of an iceberg). In the original code, a and b are both 1x25 matrices. I wanted to keep my example simple, so used 1x5 matrices instead.
K1sum = sum(K1, 'all');
Yes, I can omit the 'all' thing above and write this instead:
K1sum = sum(K1);
With the original matrices for a and b, K would be even nastier. So, yes, you are right. Using 'int' twice would take me years to compute the integral of K. So, you see I know what I am doing. The program is alright till:
K = 0.5*Ksq;
Now I need to integrate K over some fixed boundaries and I don't know how to do that with the least simulation time. Yes, I understand that singularities will exist, I will need to do something there to avoid these singularities.
"This means, IF I can assume you knoew what you were doing, before, is you just have no idea how to convery K into a function that integral2 can integrate." - Yes, you are correct here.
John D'Errico on 28 Jun 2020
You should appreciate that I can easily spend an hour or multiple hours answering just a single question on Answers. So when I see something that looks a little strange I tend to stop. What I dont want to do is spend many hours answer a question, only to have the person come back and tell me they made a mistake in line 1, and my answer now becomes irrelevant. Slow burn, acid reflux. :) You understand.
As long as what you have shown is reasonable, I would want to ask several questions (not of you) but in general.
1. What is the source of the NaNs?
2. Is the singularity integrable?
3. What is the best way to resolve the problem? Can i just isolate each singularity and treat it separately?
First, I see that each pair of elements in a and b result in one singularity. If two were close, then the singularities would be close to overlapping. That prevents me from using one solution to the problem.
a = [0, 1, 2, 3, 4];
b = [5, 6, 7, 8, 9];
Kxy(0,5)
ans =
NaN
Kxy(1,6)
ans =
NaN
As well, if we choose some random point away from the singularities, we see a relation that decreases towards zero, but as you can see, not that fast. In fact, even at the limits of your integration, the result is not that terribly close to zero.
Kxy(rand()*10,rand()*10)
ans =
0.064468
Kxy(-10,-20)
ans =
0.00036517
What have I learned so far? One option would be to split the domain into multiple subdomains. Deal with each singularity separately, possibly using some sort of Gaussian quadrature. However, these singularities are not distinct things that can easily be split apart.
What causes the NaNs? A NaN can arise from several sources. 0/0, inf-inf, inf/inf, these are just a few. We can see the cause if we look at K.
pretty(K)
/ x - 1 x - 2 x - 3 x - 4 x \2 / y - 6 y - 7 y - 8 y - 9 y - 5 \2
| ----- + ----- + ----- + ----- + -- | | ----- + ----- + ----- + ----- + ----- |
\ #4 #3 #2 #1 #5 / \ #4 #3 #2 #1 #5 /
--------------------------------------- + ------------------------------------------
2 2
where
2 2
#1 == 2 pi ((x - 4) + (y - 9) )
2 2
#2 == 2 pi ((x - 3) + (y - 8) )
2 2
#3 == 2 pi ((x - 2) + (y - 7) )
2 2
#4 == 2 pi ((x - 1) + (y - 6) )
2 2
#5 == 2 pi ((y - 5) + x )
For example, when x==1, AND y == 6, what happens? Consider what happens along some path as (x,y) approaches the point (1,6).
We have one term in there as (x-1)/#4. At some point that is within delta of the singularity, the term in question will look like
delta/(delta^2 + delta^2)
When delta is zero, this is just 0/0. For small values of delta however, it starts to have the shape delta/(2*delta^2) = 1/(2*delta). For small values of delta, that will approach infinity. Compounding this of course, is we also add to that have a second term of the form (y-6)/#4. This term also goes to infinity, and at the singularity will result in a NaN.
Next, I'll see if I can trivially catch the singularity? As we saw, near one of these singular points, Kxy fails. But a simple wrapper for this function can catch the error. I'll call it singulartrap.
singulartrap = @(X,Kmax) min(X,Kmax);
singulartrap(Kxy(1,6 + 1e-8),1e16)
ans =
1.2665e+14
singulartrap(Kxy(1,6),1e16)
ans =
1e+16
So the wrapper I created for Kxy catches both the infs and NaNs. This will allow integral2 to possibly survive.
integral2(@(x,y) singulartrap(Kxy(x,y),1e20),-A,A,-B,B)
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.
> In integral2Calc>integral2t (line 129)
In integral2Calc (line 9)
In integral2 (line 106)
ans =
14.008
>> integral2(@(x,y) singulartrap(Kxy(x,y),1e25),-A,A,-B,B)
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.
> In integral2Calc>integral2t (line 129)
In integral2Calc (line 9)
In integral2 (line 106)
ans =
16.303
>> integral2(@(x,y) singulartrap(Kxy(x,y),1e30),-A,A,-B,B)
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.
> In integral2Calc>integral2t (line 129)
In integral2Calc (line 9)
In integral2 (line 106)
ans =
22.973
>> integral2(@(x,y) singulartrap(Kxy(x,y),1e50),-A,A,-B,B)
Warning: Reached the maximum number of function evaluations (10000). The result fails the global error test.
> In integral2Calc>integral2t (line 129)
In integral2Calc (line 9)
In integral2 (line 106)
ans =
0.91402
So as I increase that maximum value, even though it is over a tiny interval, integral2 has two significant problems. First, it hits a maximum number of function evaluations, and then gets upset. But as well, changing that peak value gives us different results. And if I allow the peak to get too high, then we start to see precision problems due to the dynamic range of double precision arithmetic.
Surprisingly, integral2 does NOT offer a maximum number of function evals parameter. However, I can use other tools, for exampel quad2d, that do offer that as an option.
Warning: Reached the maximum number of function evaluations (1000000). The result fails the global error test.
ans =
22.973
While it took a lot longer to run with the high number of function evals, the result was the same as we saw before.
Regardless, this starts to beg the question of whether this function is in fact integrable at all. Is the volume under that surface a well posed problem? Or is the double integral unbounded?
I'll stop at this point, mainly because this comment has grown seriously in length, so this seems a logical break point. I'll revisit that question with another comment next.
John D'Errico on 29 Jun 2020
I'm sorry it took longer than I wanted, but sometimes it takes a fresh mind to look at a problem. But I can now show the volume under this surface is unbounded. It was not that difficult to show, but it did require me to look at the problem with a fresh mind.
In intuitive terms, you cannot compute the desired volume, because as you would work more diligently at the task, that volume essentially grows to infinity.
For a one variable problem, it becomes not much different from asking what is the integral of 1/x, from 0 to inf. The result is unbounded.
As such, I think you need to reformulate this problem, since you cannot compute any finite result.

John D'Errico on 27 Jun 2020
Edited: John D'Errico on 27 Jun 2020
Let me explain, because you either have a simple problem to solve, or an impossible one. I've restated my comment as an answer, covering either case that you may have.
While you may know what you mean by "row array" this is not something that anyone else in the world will likely know, at least not from your confusing explanation.
Is x a vector of values? So a row vector in MATLAB? It is my guess that this is what you mean. I am not at all sure if that is your intent. If it is, then you have a problem, as I would explain later if that is the case. For example, I might do this:
x = 1:5;
a = 3;
ax = x - a
ax =
-2 -1 0 1 2
In that case, ax is created as a row vector in MATLAB, possibly what you may be thinking of as a row array. We could do similarly to create a vector by. Is that what you are talking about?
Or, are these symbolic variables? For example:
syms a x b y
ax = x - a
ax =
x - a
by = y - b
by =
y - b
I am also not sure what you mean by the claim that your algebraic binomial function is not an array. I would suggest that you need to give a small example.
First, I'll consider the symbolic case, and how I would solve it trivially.
syms a x b y
ax = x - a;
by = y - b;
Now, create a binomial function of both x and y.
xyfun = ax.^2 + ax.*by - by.^3
xyfun =
(a - x)*(b - y) + (a - x)^2 + (b - y)^3
However, this is just a scalar variable. Not an array at all. Just a variable, which happens to be a function of both x and y.
I can trivially integrate over both x and y, over some arbitrary rectangular region.
int(int(xyfun,x,[0,2]),y,[0,3])
ans =
6*a^2 + 6*a*b - 21*a + 6*b^3 - 27*b^2 + 48*b - 47/2
Is that your question? If so, the answer is as easy as can be.
But, now I'll consider the other possibility, that x and y are truly each row vectors, thus "row arrays"
x = 0:5;
a = 3;
y = 0:2:10;
b = 2;
So both x and y are row vectors, I can create two new vectors, ax and by. But any relation to the original vector x is only known in tyour mind. ax is now just a vector, You need to know how that vector is related to x.
ax = x-a
ax =
-3 -2 -1 0 1 2
by = y - b
by =
-2 0 2 4 6 8
However, if you now form some new thing, from these two vectors? Each vector ax and by are really just functions of their vector index. x and y are lost. Suppose we now form the sum?
ax + by
ans =
-5 -2 1 4 7 10
This is just an other vector, treated as just a function of the vector index. At this point, any relation to the original x and y vectors is lost. You CANNOT integrate that in both x and y. It would, at this point, make no sense at all.
So either what you want to do is trivial, or it is impossible. Take your choice.

#### 1 Comment

Angshuman Sharma on 27 Jun 2020
Thank you all for trying to help me with my problem. The code that I am using is as follows:
A = 10;
B = 20;
a = [0, 1, 2, 3, 4];
b = [5, 6, 7, 8, 9];
syms x y
x = sym('x');
y = sym('y');
ax = x - a;
by = y - b;
r = ax.^2 + by.^2
K1 = - by./(2*pi.*r);
K2 = ax./(2*pi.*r);
K1sum = sum(K1, 'all');
K2sum = sum(K2, 'all');
Ksq = K1sum^2 + K2sum^2;
K = 0.5*Ksq;
syms f(x, y)
f(x, y) = @(x, y) K;
P = integral2(Ep, -A, A, -B, B)
I need to integrate K, which is a function of x and y, over their respective limits. I am so bad with coding.
When I run this code, it gives me the following error:
Error using integral2 (line 82)
First input argument must be a function handle.
Error in integration (line 53)
P = integral2(Ep, -A, A, -B, B)