3/8 Simpson's Rule

82 views (last 30 days)
Andres
Andres on 29 Sep 2022
Edited: John D'Errico on 29 Sep 2022
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
Andres
Andres on 29 Sep 2022
you should learn how to read, they are completely different questions ;)
John D'Errico
John D'Errico on 29 Sep 2022
Edited: John D'Errico on 29 Sep 2022
No. It was pretty much the same question, just asking someone to write the code for you, but providing no details of what you were trying as you did in this one. And you might be nicer to the person who just answered your question in a great deal of depth.

Sign in to comment.

Answers (1)

John D'Errico
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)
xnodes = 1×4
0 0.333333333333333 0.666666666666667 1
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)
ans =
1.71828182845905
And that is indeed exp(1) - 1, for as many digits as I recognize.
dot(fun(xnodes),[1 3 3 1])*3*h/8
ans =
1.71854015336017
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)
ans =
2
simp38(xnodes,fun)
ans =
-4.93038065763132e-32
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
xnodes = 1×4
0 3.1416 6.2832 9.4248
Between each pair of integration nodes, we now will have two new ones, equally spaced between them.
extraxnodes
extraxnodes = 1×6
1.0472 2.0944 4.1888 5.2360 7.3304 8.3776
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])
xnodes = 1×10
0 1.0472 2.0944 3.1416 4.1888 5.2360 6.2832 7.3304 8.3776 9.4248
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)
ans = 10
Effectively, we now have 3 panels for the composite simpson's 3/8 rule.
simp38(xnodes,fun)
ans = 2.0405
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)
ans = 28
format long g
simp38(xnodes,fun)
ans =
2.00038224208927
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.

Community Treasure Hunt

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

Start Hunting!