3/8 Simpson's Rule

Hello guys. I have homework and I done most of parts. But I have a problem in code. Please help me. The questions like this

Write aMATLAB implementation that applies the Simpson's 3/8 rule to find an approximation I

for I= integral(from 0 to 1) (e ^2x)*sin(6x)dx

when n=30,60,90 and 120

f=@(x)(exp^(2*x))*sin(6*x); 
 a=input('Enter lower limit a: '); 
 b=input('Enter upper limit b: ');  
 n=input('Enter the  number of sub-intervals n: '); 
 h=(b-a)/n;
if rem(n,3)~=0
   fprintf('\n Enter valid n!!!'); 
   n=input('\n Enter n as multiple of 3: ');
end
for k=1:1:n
  x(k)=a+k*h;
  y(k)=f(x(k));
end
 so=0;sm3=0;
for k=2:1:n-1
    if rem(k,3)==0
       sm3=sm3+y(k); 
     else
     so=so+y(k);
    end
end
answer=(3*h/8)*(f(a)+f(b)+3*so+2*sm3);
fprintf('\n The value of integration is %f',answer);

Answer 1

[David Hill]

f=@(x)exp(2*x).*sin(6*x); 
a=input('Enter lower limit a: '); 
b=input('Enter upper limit b: ');  
n=input('Enter the  number of sub-intervals n: ');
while rem(n,3)~=0
   fprintf('\n Enter valid n!!!'); 
   n=input('\n Enter n as multiple of 3: ');
end
x=linespace(a,b,n+1);
I=3*(b-a)/8/n*sum(f(x).*[1,3,3,repmat([2,3,3],1,(n-3)/3),1]);
fprintf('\n The value of integration is %f',I);

[Hasan Dikka]

Thank you david

Answer 2

[John D'Errico]

First, learn the basic MATLAB functions. That is, exp is used as exp(x), not exp^x. The latter makes no sense in MATLAB. exp is just a function. The problem of course, is you have gotten used to seeing things like e^x. So you assume that is how it would be written in MATLAB, or at least something vaguely like that. exp is a function, one of the most basic "special" functions you would run into, but a function that is used massively all over the place in MATLAB and in mathematics.

Next, learn to use the operators .* ./ and .^ when you will have a vector or array of elements to operate on. They will allow you to do element-wise operations. That is, if x and y are vectors, and you want the product of each element of x and y, you use x.*y.

Next, white space is important. This makes your functions more readable. There is no charge in MATLAB to add a space.

Comments are even more important. They help you to understand what you did in a spot, when you return to this code next week, next month, or next year, when all is forgotten. My suggested target for comments is there should be nearly as many lines of comments as there are code. Comments are FREE! They cost you nothing but a few seconds to type them, but they gain you much more than that when you need to debug the code. And EVERYBODY (including me) always needs to debug their code at some point in time. Next, Simpson's 3/8 rule requires 4 points per panel. It is a member of the family of Newton-Cotes rules, where we talk about a panel being a sequence of n points taken on the function, equally spaced in x.

https://en.wikipedia.org/wiki/Simpson%27s_rule

Within one panel of the rule, the weights look like

[1 3 3 1]*3*h/8

Then we can build up composite rules by creating multiple panels, one after the last. This allows us to re-use the function value from the end of the previous panel. For the Simpson 3/8 rule, the composite rule builds up the weights like this:

[1 3 3 1 + 
       1 3 3 1 +
             1 3 3 1 +
                   1 3 3 1 +
                         ...]

The final set of weights will then be:

[1 3 3 2 3 3 2 3 3 2 3 3 2 ... 3 3 1]*3*h/8

All of the closed Newton-Cotes rules work the same way, just with different sized panels, and different sets of weights. Thus we have for one panel:

Trapezoidal rule:

[1 1]*h/2

Simpson's rule:

[1 4 1]*h/3

Simpson's 3/8 rule:

[1 3 3 1]*3*h/8

Boole's rule:

[7 32 12 32 7]*2*h/45

For multiple panels, they look the same, but now the node at the joint between each panel, that node gets used twice. Therefore it effectively has double the weight.

Trapezoidal rule:

[1 2 2 2 ... 1]*h/2

Simpson's rule:

[1 4 2 4 2 4 2 4 2 ... 4 1]*h/3

Simpson's 3/8 rule:

[1 3 3 2 3 3 2 3 3 2 3 3 2 ... 3 3 1]*3*h/8

Boole's rule:

[7 32 12 32 14 32 12 32 14 32 12 32 14 ... 32 12 32 7]*2*h/45

The Simpson's 3/8 rule is just one of that family, where they all work nicely the same.

How many nodes do these rules require? Trapezoidal rule can apply to ANY number of nodes. But each panel for the basic Simpson's rule adds two more nodes. So effectively you always need an ODD number of nodes for Simpson's rule, and therefore an even number of intervals. It is usually best to think of these things in terms of panels. The Simpson's rule panel has 3 nodes in it, so it requires 2*N+1 nodes for N panels.

Similarly, Simpson's 3/8 rule uses a 4 node panel, so it requires 3*N+1 nodes, 3*N intervals, for N panels.

Only now should I try to write some code.

% Integration kernel
f = @(x) exp(2*x).*sin(6*x); 

% get the limits of integration
a = input('Enter lower limit a: '); 
b = input('Enter upper limit b: ');

% How many intervals? Anyway, rather than using input here,
% it is best to set up a loop on n. You want to have this loop over n anyway.
n = [30 60 90 120];
intapprox = zeros(size(n));
for ind = 1:numel(n)
  % ni is the number of sub-intervals. so ni/3 panels, and ni+1 nodes.
  ni = n(ind);
  
  % linspace will give ni+1 nodes
  xi = linspace(a,b,ni+1);
  
  % the stride, or node spacing is just the difference between
  % the first two elements of xi. I could also have used
  %   h = (b - a)/ni;
  % but that costs an extra divide.
  h = xi(2) - xi(1);
  
  % what is the weight vector?. We should know that ni is a multiple of 3, because
  % we are the ones who created the vector n.
  % by appending a 1 at the beginning, we start the weights properly.
  W = [1,repmat([3 3 2],[1,ni/3])];
  % The way it was done, the final weight would have been a 2. We need to reset it to 1.
  W(end) = 1;
  
  % multiply by 3*h/8 only now, after the sum is done
  intapprox(ind) = dot(f(xi),W)*3*h/8;
  
end

Did it work?

intapprox
intapprox =
     -1.0174504760826    -1.01744408724857    -1.01744374170435     -1.0174436834383
          
integral(f,0,1)
ans =
         -1.01744365644301
         
intapprox - integral(f,0,1)
ans =
  -6.81963959414666e-06  -4.30805563667036e-07   -8.5261339100029e-08   -2.6995294666321e-08

integral seems to think all is good, telling us the result is correct to about 8 decimal digits for the finer spacings.

Integral seems to have nailed the value to about as close as it could get, since we would also have used syms to perform the task.

syms X

inttrue = int(f(X),0,1)
inttrue =
3/20 - (exp(2)*(3*cos(6) - sin(6)))/20

vpa(inttrue)
ans =
-1.01744365644300762674632346283

[Hasan Dikka]

Thank you John. I solved

[Nisar]

Hi, Sir i can't make MATLAB code for new SOR-like method , please how to make it

Answer 3

[Bishal]

f=@(x)(exp^(2*x))*sin(6*x);

a=input('Enter lower limit a: ');

b=input('Enter upper limit b: ');

n=input('Enter the number of sub-intervals n: ');

h=(b-a)/n;

if rem(n,3)~=0

fprintf('\n Enter valid n!!!');

n=input('\n Enter n as multiple of 3: ');

end

for k=1:1:n

x(k)=a+k*h;

y(k)=f(x(k));

end

so=0;sm3=0;

for k=2:1:n-1

if rem(k,3)==0

sm3=sm3+y(k);

else

so=so+y(k);

end

end

answer=(3*h/8)*(f(a)+f(b)+3*so+2*sm3);

fprintf('\n The value of integration is %f',answer);