How to create loop for gaussian quadrature

clear all;

clc;

close all;

format short;

f = @(z) (200*(z/(5+z))*exp(-2*z/60))

a = 0 ;

b = 60 ;

n = 6 ;

g = input('g= ');

F = @(t) f(((b-a)*t + (b+a))/2);

His = [];

if g==2

c1 = 1 ; c2 = 1 ;

x1 = 1/sqrt(3) ; x2 = -1/sqrt(3);

G2 = ((b-a)/2)*(c1*F(x1)+ c2*F(x2));

fprintf('2 point: %.f\n',G2);

elseif g==3

c1 = 5/9 ; c2 = 8/9 ; c3 = 5/9;

x1 = -sqrt(3/5) ; x2= 0; x3 = sqrt(3/5);

G3 = ((b-a)/2)*(c1*F(x1)+ c2*F(x2)+ c3*F(x3));

fprintf('3 point: %.f\n',G3);

elseif g==4

c1 = 0.347854845 ; c2 = 0.652145155 ; c3 = 0.652145155 ; c4 = 0.347854845 ;

x1 = -0.861136312 ; x2= -0.339981044; x3 = 0.339981044; x4 = 0.861136312;

G4 = ((b-a)/2)*(c1*F(x1)+ c2*F(x2)+ c3*F(x3)+ c4*F(x4));

fprintf('4 point: %.f\n',G4);

elseif g==5

c1 = 0.236926885 ; c2 = 0.478628670 ; c3 = 0.568888889 ; c4 = 0.478628670 ; c5 = 0.236926885;

x1 = -0.906179846; x2= -0.538469310; x3 = 0; x4 = 0.538469310; x5 = 0.906179846;

G5 = ((b-a)/2)*(c1*F(x1)+ c2*F(x2)+ c3*F(x3)+ c4*F(x4)+ c5*F(x5));

fprintf('5 point: %.f\n',G5);

elseif g==6

c1 = 0.171324492 ; c2 = 0.360761573 ; c3 = 0,467913935 ; c4 = 0,467913935 ; c5 = 0.360761573; c6 = 0.171324492;

x1 = -0.932469514; x2= -0.661209386; x3 = -0.2386191860; x4 = 0.2386191860; x5 = 0.661209386; x6 = 0.932469514;

G6 = ((b-a)/2)*(c1*F(x1)+ c2*F(x2)+ c3*F(x3)+ c4*F(x4)+ c5*F(x5)+ c6*F(x6));

fprintf('6 point: %.f\n',G6);

else

disp('No solution');

end

[Jon]

What is your question?

Answer 1

[John D'Errico]

What has any of this to do with a loop anyway? There are no loops seen. You have hard coded what appear to be some basic Gauss-Legendre quadrature rules, since the nodes all seem to lie in the interval [-1,1]. (I could check that claim, but I won't bother.) DON"T WRITE CODE LIKE THAT. It will quickly lead you down the path of writing the most godawful code I can imagine.

Learn to use vectors, arrays, cell arrays, etc. You might store all of the nodes and weights in the fields of one struct, as VECTORS.

Once all of your nodes for a given order rule are stored in a vector, and the weights stored in anothe vector, then it is easy to write a loop. Or use a cell array.

Better code yet would generate the nodes and weights directly from the related orthogonal polynomials, which can be generated using recurrence relations. Or I recall a paper teaching how to use the SVD for this computation. But generating the nodes and weights on the fly is certainly beyond what you want to do. So just store the nodes and weights in some sort of array, again, as VECTORS of nodes and VECTORS of weights.

[Nguyen]

Can you make a sample of this article?

[Nguyen]

I am a beginner so please help me

[Sam Chak]

Hi @Nguyen

As a beginner, what are you expected to learn and accomplish, and how much time are you given to complete the task? It is important to understand whether you only need to acquire a portion of the coding knowledge to complete the task or if you need to borrow a MATLAB Programming book from the library.

[John D'Errico]

Store the nodes and weights like this:

GQ(1).Nodes = [0];
GQ(1).Weights = [2];

GQ(2).Nodes = [-1/sqrt(3) 1/sqrt(3)];
GQ(2).Weights = [1 1];

GQ(3).Nodes = [-sqrt(3/5) , 0 , sqrt(3/5)];
GQ(3).Weights = [5/9,8/9,5/9];

Now you can access the 1, 2, or 3 node Gauss-Legendre quadrature simply enough. You go up to the 6 node rule, but I'm not going to write it myself. That would be your job, except to point out that you don't use the full precision of a double in those coefficients. That will cost you at some point. Learn to use MATLAB.

GQ(2).Nodes

ans = 1×2

-0.5774 0.5774

GQ(2).Weights

ans = 1×2

1 1

Finally, learn to use a loop, to use the nodes and weights in that structure. Or better yet, learn to use vectorized functions, and MATLAB itself, to evaluate the function at all the nodes in one call, and then form the weighted sum of function evaluations in one line of code. Or, just use a loop. For example...

I'll integrate the function x^4+3*x^3-2*x+2, from 0 to 2, using the 1, 2, and 3 node Gaussian rules. MAKE SURE FUN IS VECTORIZED.

ftrans = @(fun,t,a,b) fun(((b-a)*t + (b+a))/2);
fun = @(x) x.^4 + 3*x.^3 - 2*x + 2;

Note that ftrans alllows the function to be passed in as an argument, as well as the lower and upper bounds of integration. So we can change those as we wish.

First, what is the exact integral? That is easy enough. We could use integral, or int.

format long g
integral(fun,0,2)

ans = 
                      18.4

Or int will work.

syms X
int(fun(X),X,0,2)

ans = 

Again, 18.4 is the answer, even though MATLAB online has written it as a fraction. Now try the Gaussian rules. An order n Gaussian rule will be exact for a polynomial of degree 2*n-1. So only the 3 node rule will be exact.

I'll write this as a helper function, to evaluate the gaussian quadrature of a function on some interval.

GQLEGINT = @(fun,order,a,b) (b - a)/2*dot(ftrans(fun,GQ(order).Nodes,a,b),GQ(order).Weights)

GQLEGINT = function_handle with value:
    @(fun,order,a,b)(b-a)/2*dot(ftrans(fun,GQ(order).Nodes,a,b),GQ(order).Weights)

GQLEGINT(fun,1,0,2)

ans = 
     8

GQLEGINT(fun,2,0,2)

ans = 
          18.2222222222222

GQLEGINT(fun,3,0,2)

ans = 
                      18.4

As predicted, that is the case. The 3 node rule was exact here. The 2 node rule was actually surprisingly good though. And the 1 node rule was, to no surprise, way out in the weeds. That will often be the case I suppose.

Finally, we can try it out on a different problem. Integrate cosine(x), from -1/2 to 1/2.

int(cos(X),X,-1/2,1/2)

ans = 

vpa(ans)

ans = 

0.95885107720840600054657587043114

And now the Gaussian rules on cos over that interval.

fun = @cos;
GQLEGINT(fun,1,-1/2,1/2)

ans = 
     1

GQLEGINT(fun,2,-1/2,1/2)

ans = 
         0.958621882624999

GQLEGINT(fun,3,-1/2,1/2)

ans = 
         0.958851569463834

Here the order 1 rule is not even that bad, but again, the order 2 and order 3 rules are much better. We got around 6 digits of precision with the order 3 rule. This is because a 4th degree polynomial quite well approximates cos(x) on that reasonably small interval. Even the 2 node rule did pretty well, since cos is pretty well approximatesd by a quadratic polynomial too.

What you should see is all of this is enabled by the ability of MATLAB to use those weights and nodes as vectors, and then a careful use of vectorization in the function fun. But really, all of this took no more than a couple of lines of code to write, once I defined the structure GQ that contains all of the nodes and weights. Now we can implement higher order rules simply by expanding GQ.

Again, if you want to use a loop, that too is easy enough. Start writing, Make an effort. If you do, you might get more help. I've shown you how to do the hard part here, because the most important part is you store the nodes and weights as vectors. All I can do i try to teach you how to use MATAB as it was designed to be used. If you really desperately think you need to use a loop, while storing all of those coefficients as scalars with numbered names, I won't help you to write bad code. So take a shot at it.

[Nguyen]

Thank you very much for helping me understand this exercise