Composite Trapezoid rule

I know for a single integral you would use the cumtrapz function. I have a problem where I have to use the composite trapz rule on a double integral and I can find nothing on how to do that. The function is

(x^2 - 2y^2 +xy^3)dxdy where the outer y limits are -1:1 the inner x limits are 0:2 and n = 2

Any help would be much appreciated. The answer is supposed to be 2.

[bym]

are you sure the answer is 2?

[Walter Roberson]

I'm not sure what n is supposed to be here. If it is the number of subdivisions to take, then perhaps the trapezoid comes out as 2 ?

[Jason]

Yes n is the number of subdivisions. The according to the book is 2 but I keep getting 8.

[Walter Roberson]

Are you remembering to multiply by the width of the interval ?

[John D'Errico]

Well, in this case, the width of the interval is exactly 1, so I'm guessing the error is in another place.

[John D'Errico]

To Jason - why not show what you have tried? You are far more likely to get help that way. I'll give you one hint. Did you start with meshgrid?

[John D'Errico]

or ndgrid?

[Jason]

This is really messy because I was forced to perform the composite trapz by hand and transfer the setup to Matlab. Doing the problem analytically and by using dblquad(for Simpson's 1/3 rule) I get 2.6667 which is the actual correct answer. Using the composite trapz rule I am supposed to get 2 with an error of 25%. Below is my attempt at the composite trapz rule for a double integral. I looked to see if there was a dblcumtrapz command but I cannot find anyway to do it in matlab, which I find kind of odd since dblquad is the Simpon's rule function command. (A) in the code below is the analytical solution which is used to find the error.

ay = -1;

by = 1;

ax = 0;

bx = 2;

n = 2;

hx = (bx -ax)/2;

hy = (by - ay)/2;

x1 = ax;

x2 = (ax +bx)/2;

x3 = bx;

y1 = ay;

y2 = (ay +by)/2;

y3 = by;

B = hy*((hx*((x1^2 - 2*y1^2 +x1*y1^3)+2*(x1^2 - 2*y2^2 +x1*y2^3)+(x1^2 - 2*y3^2 +x1*y3^3)))+2*(hx*((x2^2 - 2*y1^2 +x2*y1^3)+2*(x2^2 - 2*y2^2 +x2*y2^3)+(x2^2 - 2*y3^2 +x2*y3^3)))+(hx*((x3^2 - 2*y1^2 +x3*y1^3)+2*(x3^2 - 2*y2^2 +x3*y2^3)+(x3^2 - 2*y3^2 +x3*y3^3))))

EB = abs((A-B)/A)

[bym]

I don't see where n factors into your equations. I would suggest you look at meshgrid like John suggested

[Jason]

I'll try this again, I HAVE to you use double integration via the composite trapz rule, as for 'n' I plugged in 2 in its place for hx and hy, since h = (b-a)/n. I can not use meshgrid, and I do not see any functions for dblecumtrapz like you do with the quad function.

[Jason]

This is for a home work assignment, I have done the other 20+ problems. The prof thinks I literally have nothing else to do but sit and plug away at matlab. I have been going at this pretty solid since Thursday afternoon. The only problem I cannot get the answer to is this one.

[Jason]

For reference, this problem is part (b) of a single problem where part (a) is performed analytically with a solution of 2.6667 part (c) is solved using Simpson's 1/3 rule which can be performed with dblquad which also gives 2.6667. Part (b) has to be solved using composite trapz over the above double integral. I'm not looking for just the answer, I have the answer as 2 with a relative error of 25%. What i need to know is if there is a way to do an equivalent to a dblcumtrapz?

[bym]

why can't you use meshgrid? Is it specifically forbidden or you don't understand how to use it?

[Jason]

I guess I dont understand how using it is going to help me solve the problem using the composite trapezoid rule.

[John D'Errico]

Suppose you knew the value of this function, (x^2-2y^2+xy^3) at a grid of points. Would that help you?

[Jason]

no, not really

[bym]

how would you plot such a function?

[Jason]

I don't have the slightest idea and its not that I dont want to know, its that I dont have to plot it for this specific problem. No offense but its Sunday night, I have been doing this assignment for a few days, I have spent no time with my family,(I'm a 33 year old vet going back to school, FYI) and it really blows my mind that no where on any search engine is there any information about doing a double integral using composite trapz. Its like it has never been done before outside of this problem. At this point the few points I may lose on the assignment are not worth getting this frustrated with it.

[Jason]

If I could at least find an example I could figure it out on my own. I asked on here to see if there was some trick or function I as unaware of that might get me where I need to be. I do appreciate the help, but I have a wife to tend to and that is far more interesting than playing cat and mouse on here because no one wants to give a straight answer.

[John D'Errico]

We are not going to do your homework for you. I'm sorry that you have other things to do with your time, but if doing your homework is important to you, then you will find the time. And if it is not important to you, then why should we spend the time to do something that is not important to you anyway? There have been MANY hints offered to you.

[Jiro Doke]

To save you time, I will say that I'm not aware of any built in MATLAB function like cumtrapz for 2D. So your best bet is to implement the algorithm yourself (kind of like the way you did in one of your comments).

The command meshgrid (or ndgrid) helps you easily set up the "grid points" that you need for your algorithm. For example, notice that you went and calculated all the points based on the x and y limits and n. You could do that with meshgrid this way:

[x, y] = meshgrid(linspace(ax, bx, n+1), linspace(ay, by, n+1))

Once you have those, then it's just a matter of using those with your function "(x^2-2y^2+xy^3)" in the composite trapezoidal rule.

[Jason]

Jiro, I really appreciate that. I finally found somewhat similar application of the com trapz rule that I was able to use to figure it out. I was close originally and it was just missing a simple step. I have never used meshgrid, nor has it been discussed in class to this point so I would never have thought of using it nor would I know how to. Thank you for taking a moment to explain it. I already have the problem figured out but the meshgrid information will be useful in the future.

[Jason]

To clarify the missing part of what I had was to divide the entire B expression by 2n, which in this case would be 4. It's just too obvious, but looking at the basic form of the composite trapezoidal rule clearly shows it. Thanks again for the help even though it was something simple I missed.

Answer 1

[Sulaymon Eshkabilov]

Here is a simple solution to this exercise:

x = -1:0.02:1;
y = -1:0.025:1;
[xx,yy] = meshgrid(x,y);
Fun = xx.^2 -2*yy.^2+xx.*yy.^2;
ALL = cumtrapz(y,cumtrapz(x,Fun,2));
mesh(xx,yy,Fun)
xlabel('x')
ylabel('y')
str = '$$ \int_{-1}^{1} \int_{-1}^{1} (x^2 -2*y^2+x*y^2 )dxdy $$';
zlabel(str,'Interpreter','latex')
hold on
surf(xx,yy,ALL,'FaceAlpha',0.5,'EdgeColor','none')
plot3(xx(end),yy(end),ALL(end),'md', 'markerfacecolor', 'c')

v = [1 1 1];
view(v);