Solving Integrations with Simpson's Rule

MY CODE:

syms r
r0 = 3; % Radius of the pipe in cm

% Given equation: Q = integration of v*dA
% Cross sectional Area, A = pi*r^2
% dA = 2*pi*r*dr
% v = 2*(1-(r/r0))^(1/6)

% Declare the function
Q = 4*pi*r*(1-r/r0)^(1/6);

n = 1280; % Number of segments
a = 0; % Lower limit in min
b = 3; % Upper limit in min

simp_int = simpson(Q, n, a, b)
function integral = simpson(func, n, a, b)
     h = (b-a)/n;
     x = a;
     sum = feval(func,x);
     for i = 1:2:n-2
         x = x+h;
         sum = sum+4*feval(func,x);
         x = x+h;
         sum = sum+2*feval(func,x);
     end
     x = x+h;
     sum = sum+4*feval(func,x);
     sum = sum+feval(func,b);
     integral = (b-a)*sum/(3*n);
end

MY CODE ERROR:

Error using feval

Function to evaluate must be represented as a string scalar, character vector, or function_handle object.

Error in untitled>simpson (line 23)

sum = feval(func,x);

Error in untitled (line 19)

simp_int = simpson(Q, n, a, b)

Answer 1

[David Hill]

You are doing numerical integral, no reason for symbolics.

r0 = 3; 
f=@(r)4*pi*r.*(1-(r/r0)).^(1/6);
n = 1280; 
a = 0;
b = 3;
h = (b-a)/n;
x = a;
s=f(x);
for i = 1:2:n-2
    x = x+h;
    s = s+4*f(x);
    x = x+h;
    s = s+2*f(x);
end
x = x+h;
s = s+4*f(x);
s = s+f(b);
Integral = (b-a)*s/(3*n);

%Using built-in MATLAB functions
Integral2=integral(f,0,3);
r=0:.00001:3;
Integral3=.00001*trapz(f(r));

Answer 2

[AndresVar]

You need to declare Q as a function handle or use subs.

syms r
r0=3;
Q = 4*pi*r*(1-r/r0)^(1/6);
subs(Q,'r',0)

ans = 

0

%%% Alternatively declare Q(r)
clear;
syms Q(r)
r0=3;
Q(r)=4*pi*r*(1-r/r0)^(1/6);
Q(0)

ans = 

0

%%% feval declare Q as anonymous function
clear;
r0=3;
Q=@(r) 4*pi*r*(1-r/r0)^(1/6);
feval(Q,0)

ans = 0