Integration using function calls

How can I use the gauss quadrature function to estimate the area under a function. I also want to be able to input the number of segments used in the integration. This is what I've done so far, but I want a vector output for the gauss function that I can graph, but I keep getting just one scalar output. The first code is my main program that calls on the gauss function that is written at the bottom. My main program works well and I've tested it without the Gauss function, but something is wrong with the syntax of the Gauss function because I only get one scalar output. What's wrong?

Main Program

clear
Num = [2,4,8,10,14,18,20,30,40,50,60,80,100,1000,2500,5000];
a=-1.0;
b=12.3;
func = ('(.2*x.^3)-(2.5*x.^2)+(35.1*(sin(x).^2))+(5.23*x)-27');
for i = 1:length(Num)
  N = Num(i)+1;
      dx=(b-a)/(N-1);
      x = a:dx:b;
      f = inline(func);
      [trap_area(i)]= hw5_Trap(f,a,b,N);
      [simp_area(i)]= hw5_Simp(f,a,b,N);
      [gauss_area] = hw5_Gauss(f,a,b,N);
end
semilogx(Num,trap_area,Num,simp_area,Num,gauss_area)

Gauss Quadrature Function

function [out] = hw5_Gauss(f,a,b,N)
dx=(b-a)/(N-1);
for i = 1:N-1
   s = a + (i-2)*dx;
   e = s + dx;
   xs=((e+s)+((e-s)*(-1/sqrt(3))))/2;
   xe=((e+s)+((e-s)*(1/sqrt(3))))/2;
   gauss=(f(xs)+f(xe))*(dx/2);
end
out = gauss;

Answer 1

[Jan]

You forgot the index in the first part:

for i = 1:length(Num)
  ...
  % [gauss_area] = hw5_Gauss(f,a,b,N);
  % ==>
  [gauss_area(i)] = hw5_Gauss(f,a,b,N);
end

There seems to be something wrong with your hw5_Gauss also: The output value depends on the last iteration of the FOR loop only.

[Austin]

Right now with both of your corrections I am getting a <1x16 double> like I wanted but all the numbers are 0's and I'm pretty sure that is incorrect.

Answer 2

[Jan]

Perhaps you forgot the index:

function [out] = hw5_Gauss(f,a,b,N)
dx=(b-a)/(N-1);
out = zeros(1, N-1);
for i = 1:N-1
  s = a + (i-2)*dx;
  e = s + dx;
  xs=((e+s)+((e-s)*(-1/sqrt(3))))/2;
  xe=((e+s)+((e-s)*(1/sqrt(3))))/2;
  out(i) =(f(xs)+f(xe))*(dx/2);
end

I've used "out" directly for the output instead of writing the values at first to "gauss".

[Austin]

Yeah that did give me a vector input, but it was <1x5000 double> and I need a <1x16 double> so I can graph it on a graph with N as the N axis and the Area as the y-axis and since there are only 16 N values the vector your correction outputs is too long

[Austin]

*N as the x axix*

Answer 3

[Matt Tearle]

Your question doesn't really make sense, because the integral/quadrature approximation, for a given N, is a single number. From your reply to Jan, it seems like you're trying to plot the approximate integral as a function of N. In that case, see Jan's second answer - the problem is in the main program, not the Gauss quad. function.

However, there are some other issues. First, use function handles, not inline:

f = @(x) 0.2*x.^3 - [etc]

It's also a good idea to preallocate space for the outputs before entering the for loop:

trap_area = zeros(size(Num));

In your Gauss function, you should be doing a sum to get out! You actually don't need a loop here at all. Delete " for " and you have a vector i. Your collocation points can then be calculated in a vectorized way. Then

out = sum(f(xe) ...)*(dx/2);

[Austin]

Yeah I know I only get a single number for each given N, but I'm running the program to give a different number for all the given N's and I have 16 in the problem. Right now with both of your corrections I am getting a <1x16 double> like I wanted but all the numbers are 0's and I'm pretty sure that is incorrect.

[Matt Tearle]

Are you looking actual values, or just the plot? Because there's another issue:

Your semilog plot should be of the error, I suspect, not the actual integral values. Firstly, a semilog plot of the quadrature approximations doesn't give much insight; secondly, in this example, it gives nothing at all, because the integral is actually negative.

Now, if your quadrature routine is actually giving zeros(ie gauss_area is a 1-by-16 of zeros), then can you post your updated code, so we can see what's happening. My version of your code, with the suggestions we made incorporated, works fine.

Answer 4

[Jan]

Dear Austin, Please read our answers carefully again. The bug in hw5_Gauss has been mentioned twice now and I'm sure you have all you need to solve your homework now.

You could simplify your program noticably, e.g. you calculate "N = Num(i)+1", but access "N-1" ever. Matt's suggestion to vectorize the code would be much more efficient and matlabish. But getting the FOR-loop version to run is more important.

I suggest using the debugger: Set a breakpoint in the first line o your program and go line by line through the code to see, what happens. To learn more about the debugger, use:

docsearch debugger
docsearch breakpoint