Simpsons Rule: With for loops

Hi, So I have a question where I have to use Simpsons rule to integrate (1-x^3)*sin(x) + exp(x^2/20) between -1 and 4 with 20 intervals. The function has 4 inputs, f(x), a,b (start and end points) and n intervals

I know that I can make this code simpler with the sum function but unfortunately I have to use loops for this exercise.

My code looks like this:

function integral = simpsonsrule(f,a,b,n)
h = (b-a)/n;
x = linspace(a,b,n);
x4=0;
x2=0;
for j=2:2:b
    x4 = x4 + f(x4);
end

for k=3:2:b
    x2= x2 + f(x2);
end
integral = (h/3)*(f(a)+ f(b) + 4*(x4)+ 2*(x2));
end

And I'm calling it like this:

clear;
integral = simpsonsrule((1-x.^3)*sin(x) + exp(x.^2/20),-1,4,20)

But I'm getting the error: Undefined function or variable 'x'. but haven't I defined it with x=linspace(a,b,n)?

Answer 1

[Walter Roberson]

integral = simpsonsrule(@(x) (1-x.^3)*sin(x) + exp(x.^2/20),-1,4,20)

You need the @(x) to make an anonymous function

[Eric]

Thanks Walter.

The program works but I now see the code is wrong as I'm getting incorrect integrals. I changed my loop code to:

x4=0;
x2=0;
for j=2:2:b
    x4 = x4 + f(j);
end

for k=3:2:b
    x2 = x2 + f(k);
end

When I increase n intervals, the integral decreases. I think its calculating the area of the individual interval area but not adding them all up.

[Walter Roberson]

for j=2:2:b
    x4 = x4 + f(x(j));
end

[Eric]

No luck. Integral still keeps decreasing as n increases and doesn't stabilise. I'm thinking maybe I have to force the use of absolute values.

function integral = simpsonsrule(f,a,b,n)
h = (b-a)/n;
x = linspace(a,b,n);
x4=0;
x2=0;
for j=2:2:b
    x4 = x4 + f(x(j));
end

for k=3:2:b
    x2 = x2 + f(x(k));
end
integral = (h/3)*(f(x(a))+ f(x(b)) + 4*(x4)+ 2*(x2));
end

[Walter Roberson]

Your start and end points are numeric, not indices into x, so in your integral line you should be using f(a) and f(b)

[Eric]

Yup, changed that earlier but same result.

[Walter Roberson]

You need a small correction to the anonymous function. You currently have

@(x) (1-x.^3)*sin(x) + exp(x.^2/20)

and you need instead

@(x) (1-x.^3).*sin(x) + exp(x.^2/20)

I don't think it will change your output

[Walter Roberson]

What kind of n are you using?

[Eric]

I've been using 20,50,100,500,1000.

Answer 2

[Malek Arnous]

yes