CODEWITHSUNDEEP
matlabriemann
Cornaeleus
Cornaeleus

Reputation: 19

double integral by Riemann Sums in matlab

  clear

  n=45; // widht
  m=23; // length
  // total 990 blocks m*n


  a=-2; b=1; // x-limits
  c=2; d=4;  //  y-limits

  f=@(x,y) 4.0*x.^3.*y+0.7.*x.^2.*y+2.5.*x+0.2.*y; //function

  x=linspace(a,b,n);
  y=linspace(c,d,m);

  h1=(b-a)/n
  h2=(d-c)/m
  dA=h1*h2

  [X,Y]=meshgrid(x,y); //did a meshgrid cause q wouldnt accept different index bounds without meshgriding.
  q=sum((sum(dA*f(X,Y))))

Ive been using a formula for double riemanns at this link. https://activecalculus.org/multi/S-11-1-Double-Integrals-Rectangles.html

these are the answers

1.I=81.3000.

2.I-left=-87.4287 //-84.5523 my result

3.I-Right=-75.1072

I can't see what im doing wrong. I need input from somebody.

Upvotes: 1

Views: 1101

Answers (2)

ThomasIsCoding
ThomasIsCoding

Reputation: 101099

You didn't do anything wrong with your code. The difference comes from the resolution of x and y used in you code, since they are not sufficiently high.

For example, when you have n = 5000 and m = 5000

q = sum((sum(dA*f(X,Y)))); % your Riemann sum approach

l = integral2(f,a,b,c,d); % using built-in function for double integral to verify the value of `q`

you will see that the results are very close now

q = -81.329
l = -81.300

Upvotes: 0

Iban Cereijo
Iban Cereijo

Reputation: 1675

I would debug the integration scheme with a dummy constant function 

f = @(x,y) 1.0*ones(size(x))

The result should be the exact total area (b-a)*(d-c) = 6, but your code gives 6.415.

The problem is that you are integrating outside the domain with those limits. You should stop the domain discretization one step before in each dimension:

h1 = (b-a)/n
h2 = (d-c)/m
x = linspace(a, b - h1, n);
y = linspace(c, d - h2, m);

This will give you the expected area for the dummy function:

q =  6.0000

and for the real function, evaluating at the top-left corner, you get:

q = -87.482

Upvotes: 0

Related Questions