# 3/8 Simpson's Rule

35 views (last 30 days)

Show older comments

Hello, i have to solve an integral ((x^2+exp(2*x))*(sin(x))) but the code must be iterative, it means that, for each iteration the value of the integral must get closer to the real value. I also have to calculate the absolute and relative error, but when executing the code, the table doesn't calculate the value for each interval, it just shows the same answer for all of them.

clear all;

close all;

clc;

format long

fprintf('REGLA DE SIMPSON 3/8 \n');

fprintf('Este programa usa la regla de Simpson 3/8 para aproximar el valor de una integral definida \n');

syms x

f=(x^2+exp(2*x))*(sin(x)); % Función

a=-1; % Intervalo inferior

b=1; % Intervalo superior

n=input('Ingrese el número de iteraciones: ');

fprintf('Iteración Solución Error real Error aprox\n')

for k=1:n

while(mod(n,3)~=0)

n=n+1;

end

h=(b-a)/n;

suma=subs(f,a)+subs(f,b);

for i=1:n-1

if(mod(i,3)==0)

suma=suma+2*subs(f,a+i*h);

valor_experimental=double(suma*3*h/8);

valor_teorico=double(int(f,a,b));

error_abs=(valor_experimental-valor_teorico);

error_r=double(abs(((n*h^5)/80)*max(subs(diff(f,4),a),subs(diff(f,4),b))));

else

suma=suma+3*subs(f,a+i*h);

valor_experimental=double(suma*3*h/8);

valor_teorico=double(int(f,a,b));

error_abs=(valor_experimental-valor_teorico);

error_r=double(abs(((n*h^5)/80)*max(subs(diff(f,4),a),subs(diff(f,4),b))));

end

end

fprintf('%6.2f %6.5f %0.5e %0.5e\n',k,valor_experimental,error_abs,error_r);

end

fprintf('El valor teórico de la integral es: %6.5f',valor_teorico);

##### 8 Comments

John D'Errico
on 29 Sep 2022

Edited: John D'Errico
on 29 Sep 2022

### Answers (1)

John D'Errico
on 29 Sep 2022

Edited: John D'Errico
on 29 Sep 2022

Please don't keep on posting the same question. You got no answer, since your question was itself terribly confusing, as was your code. Sorry. It was. However, you did make some fairly credible effort, so I'll show how you might solve it on a related problem. First, what is Simpson's 3/8 rule?

In there, we see it is a higher order Newton-Cotes rule. Well, higher order than the standard Simpson's rule. The idea is we have "panels" of 4 nodes, with weights that look like

(f(0) + 3*f(h) + 3*f(2*h) + f(3*h))*3*h/8

So a [1 3 3 1] pattern, with 3h/8 as a scale factor on the outside.

It approximates the integral of a function over the interval [0,3*h]. And of course, you can have multiple panels side by side. How well does it work? For example:

fun = @exp;

a = 0;

b = 1;

h = (b-a)/3;

format long g

xnodes = linspace(0,1,4)

So that is ONE panel. (Panel is the standard name used here.) Lets try it out. First, we know the correct integral over that domain is just exp(1)-1.

integral(fun,0,1)

And that is indeed exp(1) - 1, for as many digits as I recognize.

dot(fun(xnodes),[1 3 3 1])*3*h/8

which is indeed pretty decent, with about 4 significant digits correct. Now what does an iterative version of this rule look like? For this, we first need to build a Simpson's 3/8 tool, that will work on multiple panels side by side. Note that the node shared by each panel gets twice the weight, just as in Simpson's rule or trapezoidal rule. So I'll pass in a set of nodes. No need to pass in h also, as that can be found from the nodes themselves. I'll also pass in the function to be evaluated. Note that this code is inefficient, because it will not reuse the function values from a previous iteration. Life sucks. Better code would be efficient, but it would also be more confusing to a novice.

I'll try it now. First, I'll make the problem slightly harder, computing the integral of sin(x) over 1.5 periods, so [0,3*pi]. We expect this integral to be 2. (I've made it a slightly harder problem, because I want you to see the convergence after a couple of iterations.)

xnodes = linspace(0,3*pi,4);

fun = @sin;

integral(fun,0,3*pi)

simp38(xnodes,fun)

This time, I would expect to see complete crap for a result. In fact, essentially zero is what I would expect for a one panel rule here, but h is pretty darn large, at pi/2 for this first case.

But now is when we start to think about how to make the solution iterative. First, we have xnodes. Next, between each of those nodes, I'll insert two NEW nodes. It will look like this:

extraxnodes = sort([xnodes(1:end - 1)*2/3 + xnodes(2:end)*1/3 , xnodes(1:end - 1)*1/3 + xnodes(2:end)*2/3]);

Now the first panel had nodes where?

format short

xnodes

Between each pair of integration nodes, we now will have two new ones, equally spaced between them.

extraxnodes

We can plot the nodes, using different heights so you can see what to expect.

plot([1;1]*xnodes,repmat([0;1],1,4),'b-o',[1;1]*extraxnodes,repmat([0;1/2],1,6),'r-s')

We can visualize the first panel as the blue nodes. Then I inserted 2 points between each original pair of nodes. So now we have 3 panels, and 10 total points. Between every pair of points in the original rule, I now have a complete new panel. You should see the nodes between each panel are essentially shared. This is a standard thing with all newton-cotes composite rules. (Look back at the trapezoidal rule, or the composite standard Simpson's rule. You will see the same fundamental idea.)

So the new rule now has

xnodes = sort([xnodes,extraxnodes])

I could have done that more elegantly I guess, but again, life sucks. Do you see that now, we have a new set of nodes, with the new nodal spacing exactly 1/3 of the previous one? Now instead of 4 nodes, we should expect to see 10 nodes.

numel(xnodes)

Effectively, we now have 3 panels for the composite simpson's 3/8 rule.

simp38(xnodes,fun)

As you can see, now we are coming reasonably close to the target of exactly 2, as we know the true integral must be. Were I to repeat the above process, inserting 2 more nodes between each pair of old nodes, we would then have 9 total panels, and 28 points overall, if I did my caclulations correctly.

extraxnodes = sort([xnodes(1:end - 1)*2/3 + xnodes(2:end)*1/3 , xnodes(1:end - 1)*1/3 + xnodes(2:end)*2/3]);

xnodes = sort([xnodes,extraxnodes]);

numel(xnodes)

format long g

simp38(xnodes,fun)

As you can see, the solution is converging on 2. With some thought, you could derive the convergence rate, as h decreases. This third iteration has a nodal spacing of pi/18, so we went from a nodal spacing of pi/2, to pi/6, to pi/18.

function predint = simp38(xnodes,fun)

% compute Simpson's 3/8 estimate of the integral over the support of xnodes

% get h

h = xnodes(2) - xnodes(1);

n = numel(xnodes);

% n MUST be an integer of the form 3*Panels + 1, else the composite

% Simpson's 3/8 rule must fail.

if (n < 4) || (mod(n,3) ~= 1)

error('Help! The sky is falling! Wrong number of nodes in xnodes.')

end

% the weights are simple. They will look like 1 3 3 2 3 3 2 3 3 2 ... 2 3 3 1

npanels = (n-1)/3;

W = repmat(3,1,n);

if npanels > 1

W(4:3:end) = 2;

end

W([1,end]) = 1;

W = W*3*h/8;

% evaluate the function at all nodes. (Make sure fun is vectorized!)

ynodes = fun(xnodes);

% form the composite Simpson's 3/8 rule estimate.

predint = dot(ynodes,W);

end

It would be my guess that this is your target, to write code that would iteratively subdivide the intervals, taking care to do so to be consistent with a Simpson's 3/8 rule. Note that I have not solved your explicit problem, as that would be doing your homework. As it is, I have come perilously close to doing that for you here. But you will need to do some work to solve your specific problem, even if you used my code. You probably want to write a loop, that will iteratively refine the intervals.

And of course, if you did use my direct code, your teacher might wonder just a bit, since it would be a bit more sophisticated for what they would expect from a student. So you need to do some thinking, make some extra effort. Such is life. In a just world, the rewards are often correlated with the effort expended. Regardless, I would strongly recommend reading this answer at least twice. Think about what each line of code I wrote does, and why it works, but most importantly, what it does in the end.

##### 0 Comments

### See Also

### Categories

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!