gaussian quadrature

[Antonio]

Hello,

First of all sorry for my bad English .
I'm making a program to calculate integrals by gaussian quadrature and I have the following problem , I've done two different ways.
Program 1 , due to the unique character of this calculator , sometimes it works, sometimes not , so i try to make the program 2
Program 2 , do not put the correct result in a0 .

Program 1

Code:

#pragma mode( separator(.,;) integer(h64) ) 


EXPORT Cuadratura_Gaussiana()
// Calculo de las raices de los polinomios de Legendre
BEGIN
    local i,aux;
    local nr:=3;
    M1:=makemat(0,nr,2); // 1ª col factores de ponderación ,2ª col raices
    // Cálculo de las raices  de los polinomios de Legendre
    M2:=CAS.zeros(legendre(nr)); 
    
    // Cálculo de los coeficientes 
    F1:=CAS.diff(legendre(nr));

    for i from 1 to nr do
         M1(i,2):=M2(i);
         aux:=F1(M2(i));
         M1(i,2):=2/((1-M2(i)^2)*(aux^2)); // coeficiente
    end;
     
end;

Programa 2

Code:

#pragma mode( separator(.,;) integer(h64) )

coeficientes();

EXPORT Cuadratura_Gaussiana()
BEGIN
 
    local n:=3;
    M1:=makemat(0,n,4); 
    
    coeficientes(n);    
    msgbox("Orden de las columnas: c_i, r_i , F(r_i), c_i*F(ri)");
    editmat(M1);
end;

#cas
coeficientes(nr):=
BEGIN
    // Cálculo de las raices y los coeficientes para la cuadratura gaussiana de de nr términos
    local i,a0;
    M2:=zeros(legendre(nr)); //raices del polinomio
    // Cálculo de los coeficientes de ponderacion
    for i from 1 to nr do
         M1(i,2):=M2(i);//  raiz del polinomio; 
         a0:=subst(diff(legendre(nr)),x=M2(i));
         L1(i):=a0; 
         M1(i,1):=2/((1-M2(i)^2)*(a0^2));//coeficiente;
    end;
    
END;
#end

Let's see if someone can lend me a hand, thanks .

[Didier Lachieze]

(08-23-2015 07:01 PM)Antonio Wrote:  Program 1 , due to the unique character of this calculator , sometimes it works, sometimes not , so i try to make the program 2

Can you explain what happens when it doesn’t work? On my side it seems to works fine.

(08-23-2015 07:01 PM)Antonio Wrote:  Program 2 , do not put the correct result in a0 .

I’ve done some tests with your program 2 and if I’m correct it seems that there is an issue when the CAS function legendre() is called from a CAS program:

• When called from CAS it returns a polynom of x, for ex. legendre(3) returns 5/2*x^3-3/2*x
• However when called from a CAS program it returns a polynom of X instead of x: 5/2*X^3-3/2*X
  
  For example if you add the following line in your code: f:=legendre(nr); just before a0:=subst(diff(legendre(nr)),x=M2(i)); then after running your program, if you type CAS.f you will get 5/2*X^3-3/2*X.
  Same for diff(legendre(nr)) which returns: 15/2*X^2-3/2
  As X is a global predefined variable, then a0 is always evaluated to the same value and just depends on the value of X (-3/2 if X=0, its default value).

[Antonio]

Hello Didier,

Thanks for the reply.

Regarding the program 1 , what I mean is that sometimes it runs and calculates the roots and coefficients and sometimes fails to execute
M2:=CAS.zeros(legendre(nr)); or
F1:=CAS.diff(legendre(nr));
does not give any error message unless it does not.

Regarding the program 2, I tried to change x by X , as you suggested, I had to add one more line to assign the value to the X ,as it does not properly
a0:=subst(diff(legendre(nr)),X=M2(i));

they work as

X:=M2(i);
a0:=subst(diff(legendre(nr)),X=0);

Another problem is that not running the loop, and I 've been manually deactivated by assigning values ​​to the variable i, work when i = 2 or i = 3 , but not when i = 1, the error is in

M1(i,1):=2/((1-M2(i)^2)*(a0^2));
and I dont see the failure more I look.

[Didier Lachieze]

(08-25-2015 09:23 AM)Antonio Wrote:  Regarding the program 1 , what I mean is that sometimes it runs and calculates the roots and coefficients and sometimes fails to execute
M2:=CAS.zeros(legendre(nr)); or
F1:=CAS.diff(legendre(nr));
does not give any error message unless it does not.

On my side every time I've tried the program 1 it has worked correctly so I can't reproduce your issue.

(08-25-2015 09:23 AM)Antonio Wrote:  Another problem is that not running the loop, and I 've been manually deactivated by assigning values ​​to the variable i, work when i = 2 or i = 3 , but not when i = 1, the error is in

M1(i,1):=2/((1-M2(i)^2)*(a0^2));
and I dont see the failure more I look.

I think this is because M1 is predefined as a matrix of reals and you're trying to assign to it an exact number when a0 is -3/2. Try with adding a dot to the first 2 to get an approximate number (or create your own matrix variable, different from the predefined ones M0-M9):

M1(i,1):=2./((1-M2(i)^2)*(a0^2));

[Antonio]

Already solved

As you said I put the point to take the number as Real rather than exact and have solved two problems , which did not work when i = 1 and the problem of the loop. Although I do not see the relationship between the loop and the latter problem that the variable a0 can be real or exact.


The corrections program that works is

Code:

#pragma mode( separator(.,;) integer(h64) )

coeficientes();

EXPORT Cuadratura_Gaussiana()
BEGIN
 
    local n:=3;
    M1:=makemat(0,n,4); 
    
    coeficientes(n);    
    msgbox("Orden de las columnas: c_i, r_i , F(r_i), c_i*F(ri)");
    editmat(M1);
end;

#cas
coeficientes(nr):=
BEGIN
    // Cálculo de las raices y los coeficientes para la cuadratura gaussiana de de nr términos
    local i,a0;
    M2:=zeros(legendre(nr)); //raices del polinomio
    // Cálculo de los coeficientes de ponderacion
    for i from 1 to nr do
         M1(i,2):=M2(i);//  raiz del polinomio;
         X:=M2(i);
         a0:=subst(diff(legendre(nr)),X=0);  
         M1(i,1):=2.[/b]/((1-M2(i)^2)*(a0^2));//coeficiente;
    end;
    
END;
#end

In relation to program 1 , I have not noticed any particular pattern in failure . It happens also in other programs which use CAS.funtion( ) sometimes makes the call to CAS successfully and some not.

Every time I like least this calculator , since to do something that comes out of the most basic , even as simple as this program is a headache..

In any case thank you very much for your help Didier