Why does my function give weird outputs??

I am trying to create a code to calculate work. I have a vector for n and when I run the script I keep getting numbers that don't make sense such as a zero for the estimated integral when using the vector. when I run the script without the n vector and just input n I get more logical answers but they still seem off.

function [ w,Ws,error ] = trap_wfp610( ~ )
% Calculates the work done using numerical integration and the trapezoidal rule
% Detailed explanation goes here
n=[1,11,21,3,1,41,51]';
h=(5*10^-2)/n;
a=0;
b=5*10^-2;
f=@(a)(5*(10^5)*a+10^5)*pi*(5*10^-2)^2;
g=@(b)(5*(10^5)*b+10^5)*pi*(5*10^-2)^2;
w=h/2*(f(a)+g(b));

syms x
k1=500000;
k2=10^5;
p=(k1*x+k2)*pi*.05^2;
Ws=int(p,x,0,5*10^-2);

error=abs((w-Ws)/Ws)*100;

disp('Work Trap Rule');
disp(w);
disp('Work Actual');
disp(Ws);
disp('Error');
disp(error);

plot(w,n)
plot(error,n)
end

These are my outputs

>> trap_wfp610
Work Trap Rule
         0         0         0         0         0         0    0.8662

Work Actual

(225*pi)/16

Error

[ 100, 100, 100, 100, 100, 100, 5000/51]

ans =

           0         0         0         0         0         0    0.8662

which seem illogical.

[John D'Errico]

Please learn to format your code so it becomes readable.

Answer 1

[John D'Errico]

Your post gives weird outputs because you do not understand many basics of MATLAB, from what I can see so far. For example...

n=[1,11,21,31,41,51]';
h=(5*10^-2)/n;

So you create a vector, then do a divide but you need to check that what you thought you did and what you actually did were two totally different things. You need to learn how to use element-wise division, not just use the / operator. That is, use ./ instead. There is a difference.

So, lets see what happens when you execute the first two lines of your code.

h
h =
     0    0    0    0    0   0.00098039

Not what you expected, I suppose. The tutorials teach that element-wise division (and often multiplication) requires the use of the ./ operator.

h=(5*10^-2)./n
h =
         0.05
    0.0045455
     0.002381
    0.0016129
    0.0012195
   0.00098039

Forward slash does a linear algebra operation, that is not what you expected to see. The way I did it does what you expected to see.

I suggest you learn what the basic operators do, test each line of code to verify that you got what you expected.

Next, it looks like you are trying to do an n-point trapezoidal rule, with n varying over the set [1,11,21,3,1,41,51]. A good idea on the face of it, but you didn't generate any set of points on that interval, no loops, no implicit loops in your code. You just subtracted two numbers. MATLAB does not know what you are trying to do here. How is MATLAB supposed to interpret n?

There are other things I could comment on. For example, while I see while you have learned the basic syntax of function handles, my guess is you do not understand them. In your code, you created f and g as identically the same functions, one is a function of a, and one a function of b. But inside, they are identical.

I suggest you start learning the basics. For example, suppose you want to perform a basic trapezoidal rule on the function x.^2, over the interval [0,1]. See that I used .^, to allow the operation to work for a vector of inputs, an element-wise exponentiation. (Addition and subtraction do not need element-wise operators. * and / also sometimes do what you will have expected, but until you learn how they work, it is best to be careful.)

f = @(x) x.^2;
n = 11;
xi = linspace(0,1,n);
yi = f(xi);

trapz(xi,yi)
ans =
                   0.335

Not too bad, since 1 / 3 is the correct answer. The 11 point trapezoidal approximation seems about right. You could do the trapezoidal rule by computing the weights, and writing it as a weighted sum

Be careful though, since I'll bet that you think you could have used

n = [1 11 21 31 41 51];

in the above example, and that MATLAB would have known you wanted to do the same operation for each value of n. Sorry, but MATLAB will get that wrong.

I could write more sophisticated code that would compute n point trapezoidal rule estimates for each of those values of n, but it would just introduce a nifty trick that might be confusing.

You need to spend some time learning how to use MATLAB. Learn the basics first.

[Pedram Safaee]

Very arrogant comment. This should be flagged for removal so people do not think that they can bulley others just because they supposedly know something better than others.