Problem in generating iterations in 3/8 simpson rule

I want to make simpson rule algorithm for my function for different intervals (here, 3^p), and plot the error vs p. The code below is showing error in the use of syms in 1st line and the iteration of line 9 (Itr=Itr+f(N+1)/2). Please tell what is wrong in this and how can I correct it.

syms p
N=3^p;
h=1/N;
    for i=1,N;
        x=i*h;
        f(i)=exp(-x^2)
        Itr=f(1)/2;
        Itr=Itr+(3/2)*f(i);
        Itr=Itr+f(N+1)/2;
    end
    Itr=Itr*h;
for p=1,10;
plot(p,itr)
end 

Error using sym/subsindex (line 857)
Invalid indexing or function definition. Indexing must follow MATLAB indexing. Function arguments must be symbolic variables, and function body must be sym expression.

Error in sym/subsref (line 902)
            R_tilde = builtin('subsref',L_tilde,Idx);

Error in simpson3by8rule (line 9)
        Itr=Itr+f(N+1)/2;

Answer 1

[Alan Stevens]

Not entirely clear to me between what limits you are integrating, nor what your loops are doing! Perhaps the following might help:

% Integrate exp(-x^2) from x = 0 to x = 3 using Simpson's 3/8 rule 
% I = 3*h/8*(f(a))+sum(f(i)+f(i+1),i=2:3:n-2)+sum(f(i),i=3:3:n-1)+f(b))

f = @(x) exp(-x.^2);   % function
a = 0;                 % lower limit
b = 3;                 % upper limit

n = 333;         % number of panels (multiple of 3)
h = (b - a)/n;   % panel size
x = a:h:b;       % vector of x values
  
s3 = 0;         
for j = 2:3:n-2
    s3 = s3 + f(x(j)) + f(x(j+1)); 
end
s2 = 0;
for j = 3:3:n-1
   s2 = s2 + f(x(j)); 
end

I = 3*h/8*(f(a) + 3*s3 + 2*s2 + f(b));
disp('Approximate integral of exp(-x^2) from x=0 to x=3 using')
disp(['Simpson''s 3/8 rule with ' int2str(n) ' panels is: '  num2str(I,4)])

Modify to suit exactly what you want.

[Ipshita]

Sorry for not specifying my question properly. I have to integrate the given function from 0 to 1 in 3^p intervals where p varies from 1 to 10. Then I have to plot graph between the error and p. So the doubt is that how can I represent p here?

[Alan Stevens]

Ony a small modiication required, using p as a loop index:

Itrue = 0.7368;
f = @(x) exp(-x.^2);   % function
a = 0;                 % lower limit
b = 1;                 % upper limit

I = zeros(10,1);       % preallocate space
for p = 1:10           % Loop ten times
    
    n = 3^p;         % number of panels (multiple of 3)
    h = (b - a)/n;   % panel size
    x = a:h:b;       % vector of x values

    s3 = 0;          
    for j = 2:3:n-2
        s3 = s3 + f(x(j)) + f(x(j+1)); 
    end
    s2 = 0;
    for j = 3:3:n-1
       s2 = s2 + f(x(j)); 
    end

    I(p) = 3*h/8*(f(a) + 3*s3 + 2*s2 + f(b));
    
end

errors = I - Itrue;

plot(1:p,errors,'-o'),grid
xlabel('p'),ylabel('error')

[Ipshita]

Thankyou very much for your answer. As I am new to MATLAB, can you please explain the following:

1.plot(1:p,errors,'-o')

Why have you taken 1:p here, and not just p? When I take plot(p,errors,'-o'), why is the graph plotted just for the value p=10?

2. f = @(x) exp(-x.^2)

@ is used here to denote function handle? I read about it but did not understand exactly. Without using @(x), the command window is showing this:

Array indices must be positive integers or logical values.

Error in simpson3by8 (line 11)

s3 = s3 + f(x(j)) + f(x(j+1));

Thanks in advance!

[Alan Stevens]

p is a counter, ending up with the single value of 10. Hence, we need to use a vector 1:10 (or 1:p) as the "x" values in the plot command.

Yes, @ creates a pointer to a function.

[Ipshita]

Thankyou!

[Ipshita]

Hello, I have to plot both the error as well as the estimate error using a formula. But, it is showing an error as below. Can you please tell what is the mistake and how to correct it? Below is just the starting part, and not the whole code.

f = @(x) exp(-x.^2);  
a = 0;                 
b = 1;                 
for p=1:20
N=3*p; % no. of intervals
h=1/N;
x = a:h:b;       
df4 = matlabFunction(diff(sym(f), 4))
estimateerror(x)=-3*df4(x)*h^5/80

df4 =

  function_handle with value:

    @(x)exp(-x.^2).*1.2e+1-x.^2.*exp(-x.^2).*4.8e+1+x.^4.*exp(-x.^2).*1.6e+1

Array indices must be positive integers or logical values.

Error in sympson3by8attempt4 (line 9)
estimateerror(x)=-3*df4(x)*h^5/80

[Alan Stevens]

In

estimateerror(x)=-3*df4(x)*h^5/80

the first value of x is 0, but indices in Matlab start at 1 not 0. Try replacing it with

for i = 1:numel(x);
    estimateerror(i)=-3*df4(x(i))*h^5/80;
end

(Matlab would also not like other values of x that were not integers being used as indices).

[Ipshita]

Thanks a lot, it worked!