I am trying to do a double integral of a function using the simpsons 1/3 rule, need help with setup

I am just needing to know how to set the variables up. This is the code i have so far,

n = 2;
func = @(x,y) x^3 + 2*x*y-y^2; 
x = a;
h = (b-a)/n;
m = n/2;
s = 0;
s2 = 0;
for i=1:m
  s = s + func(x) + 4*func(x+h) + func(x+2*h);
  x = x + 2*h;
  s2 = s2 + func(y) + 4*func(y+h) + func(y+2*h);
  y = y + 2*h;
end
I = (b-a)*s/(3*n);
I = (b-a)*s2/(3*n);

[EDITED, Jan, code formatted to make it readable]

[Cam Salzberger]

Hello Garnet,

I'm a little confused by your code. You assign n = 2, m = n/2, and then loop from 1 to m (when m is 1). Why do you have a loop? Also, the last two lines both assign values to "I". Did you mean to assign to "I2" on the last line, or something like that?

Answer 1

[mohammed almashri]

function I =simsonss(func,a,b,n)
x=a; h=(b-a)/n;
s=func(a);
for i=1:2:n-1
    x=a+i*h;
    s=s+(4* func(x));
end
for j=2:2:n-2
    x=a+j*h;
    s=s+2*func(x);
end
s=s+func(b);
s=s*(b-a)/(3*n);
end

Answer 2

[Zoe]

this is an example for calculating double integrals with simpsons rule, pretty simple I think:

f = @(x,y) x+y.^3;
a = 0;
b = 4;
c = 1;
d = 3;
n = 4;
m = 2;
s = doubleint(f,a,b,c,d,n,m);

s = 0

s_ana = 0.5*(b^2-a^2)*(d-c)+0.25*(b-a)*(d^4-c^4)

s_ana = 96

disp('the approximate integral of Simpson rule is:')

the approximate integral of Simpson rule is:

disp(s)

    96

disp('analytically:')

analytically:

disp(s_ana)

    96

function s = doubleint(f,a,b,c,d,n,m)
hx = (b-a)/n;
hy = (d-c)/m;
s = 0
for i = 0:n
    if i == 0 || i == n
        p = 1;
    elseif rem(i,2) ~= 0 
        p = 4;
    else 
        p = 2;
    end
    for j = 0:m
        if j == 0 || j == m
            q = 1;
        elseif rem(j,2)~= 0 
            q = 4;
        else 
            q = 2;
        end
        x = a + i*hx;
        y = c + j*hy;
        s = s + p*q*f(x,y);
    end
end
s = (hx*hy)/9 * s;
%disp('the approximate integral of Simpson rule is:')
end