# Solving Integrations with Simpson's Rule

5 views (last 30 days)
Mckale Grant on 23 Feb 2022
Edited: AndresVar on 23 Feb 2022
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)

David Hill on 23 Feb 2022
Edited: David Hill on 23 Feb 2022
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));

AndresVar on 23 Feb 2022
Edited: AndresVar on 23 Feb 2022
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