Index must be positive integer or logical. I keep getting the same error when trying to run my code for numerical integration

3 views (last 30 days)
Hellooo, I keep getting the same error in my code for numerical integration. The code has a default program and also user inputs where the user can put in values for the limits and a function and then to plot the integration rules on a graph. I keep getting the same error 'Attempted to access fx(0); index must be a positive integer or logical.' I think the problem is with what i'm defining 'fx' as because it is a vector instead of a function and I really can't think of a way to solve the problem. I have tried using str2func and the 's' input as well. Help with this problem or if someone could show me what i'm doing wrong that would be fantastic! My code is:
% Implentation of Rectangle, Trapezium and Simpson's Rule % a: lower bound, b: upper bound, N: number of strips
clear, %Default program for run = menu('Run default program?','Yes','No'); if run == 1; a = 0; b = 10; N = 50; h = (b-a)/N; x = linspace(a,b); fx = (exp(-x/4)).*cos(x); elseif run == 2; fprintf(' \n Numerical Integration'); fprintf(' \n Input values for integral boundaries a and b with a<b') fprintf(' \n a = '); a= input(''); fprintf(' \n b = '); b= input('');
if a > b;
fprintf(' \n a has to be less than the value of b')
fprintf(' \n Input new values for a and b')
fprintf(' \n a = '); a= input('');
fprintf(' \n b = '); b= input('');
end
fprintf(' \n Input the number of strips')
fprintf(' \n N = '); N= input('');
if mod(N,2)> 0;
N = N+1;
else
end
x = linspace(a,b);
f = input('type in f(x) = ');
fx = @(x) feval(f, x);
h = (b-a)/N;
end
end
% Rectangle Rule
rectsum = 0;
for i = 1:N-1
rectsum = rectsum + fx(i);
end
Rectangle_Rule = h*rectsum;
% Trapezium Rule
trapsum = 0;
for i = 1:N-1;
trapsum = trapsum + 2*fx(i);
end
Trapezium_Rule = h/2*(fx(a)+ trapsum + fx(b));
% Simpson's Rule
simp4 = 0;
simp2 = 0;
for i = 1:2:N-1
simp4 = simp4 + 4*fx(i);
end
for i = 2:2:N-2
simp2 = simp2 + 2*fx(i);
end
Simpsons_Rule = h/3*(fx(a)+ simp4 + simp2 + fx(b));
% Exact Value
pp = spline(x,fx);
Exact_Value = integral(@(x)ppval(pp,x), a, b);
% Plotting the graph
clf, hold on;
plot(x,fx);
xlabel('x', 'FontWeight', 'bold');
ylabel('fx', 'FontWeight', 'bold');
title('Numerical integration of function fx' , 'FontWeight', 'bold');
xstring = num2str(Exact_Value);
text(1.5,20,['Exact Value =' xstring]);
xstring2 = num2str(Simpsons_Rule);
text(1.5,18,['Simpsons Rule =' xstring2]);
xstring3 = num2str(Trapezium_Rule);
text(1.5,16,['Trapezium Rule =' xstring3]);
xstring4 = num2str(Rectangle_Rule);
text(1.5,14,['Rectangle Rule =' xstring4]);
axis 'auto xy'
% Plot Rectangle Rule
hold on
for x=a:h:(b-h);
left = x; right = x+h; bottom = 0; top = fx(x);
X = [left left right right]; Y = [bottom top top bottom];
fill(X,Y, 'b', 'FaceAlpha', 0.3)
end
% 9. Display a smooth line in the graph
x = a:h/100:b;
plot(x, fx(x), 'k')
  1 Comment
Justin
Justin on 29 Apr 2014
Hi
I wasn't able to get a good copy of the code to paste to test. I don't want to experiment with it and guess where the if/end statements are correct.
If you can try to edit the question to format it better that may help. Make sure you look at the preview of your text before submitting it to see if it looks right.

Sign in to comment.

Answers (1)

Geoff Hayes
Geoff Hayes on 29 Apr 2014
Hi Beth,
The error message
'Attempted to access fx(0); index must be a positive integer or logical.'
arises because in MATLAB the indices for an array/vector of length N are 1,2,…,*N*. So the first element is always accessed with a 1 i.e. fx(1).
In your default program case, you set a to zero and then you (later) try to evaluate f(a) for the Trapezium Rule and Simpson's Rule. Remember that your vector fx contains the results of evaluating some function in the interval [a,b] so you want to access the result of evaluating the function at a which is fx(1). The equivalent of fx(a) and fx(b) would be fx(1) and fx(end) respectively.
Geoff
  2 Comments
beth
beth on 29 Apr 2014
Ah thank you Geoff! I managed to process the graph and it looks good however the values I get for the Simpsons, Rectangle, Trapezium and Exact are all wrong. I was wondering if you knew why this was the case when I changed fx(a) to fx(1) and and fx(b) to fx(end)?
Geoff Hayes
Geoff Hayes on 29 Apr 2014
Looking at the equation outlined in Trapezium Rule what you have written seems almost correct. Your code for the Trapezium Rule is probably now something like:
Trapezium_Rule = h/2*(fx(1)+ trapsum + fx(end));
Your h=(b-a)/N and then the division by two seems fine, but I think that you should review trapsum. Remember that MATLAB arrays start at one (not zero) so I think the code is adding the first element fx(1) too many times. Just start from i=2 and try again. You may want to review for the other equations too.

Sign in to comment.

Categories

Find more on Search Path in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!