Recent

Author Topic: Gaussian Quadratures  (Read 1664 times)

user07022024

  • New Member
  • *
  • Posts: 13
Gaussian Quadratures
« on: July 12, 2024, 01:54:35 pm »

I found that following procedure from numerical recipes can be useful for calculating nodes and weights

Code: [Select]

TYPE
   glnp = ARRAY  OF real;
   glnpnp = ARRAY  OF glnp;


FUNCTION sign(a,b: real): real;
   BEGIN
      IF (b < 0) THEN sign := -abs(a) ELSE sign := abs(a)
   END;

PROCEDURE tqli(VAR d,e: glnp; n: integer; VAR z: glnpnp);
(* Programs using routine TQLI must define the types
where np is the physical dimension of the matrix to be analyzed. *)
LABEL 1,2;
VAR
   m,l,iter,i,k: integer;
   s,r,p,g,f,dd,c,b: real;
BEGIN
      FOR i := 2 TO n DO BEGIN
         e[i-1] := e[i]
      END;
      e[n] := 0.0;
      FOR l := 1 TO n DO BEGIN
         iter := 0;
1:         FOR m := l TO n-1 DO BEGIN
            dd := abs(d[m])+abs(d[m+1]);
            IF  (abs(e[m])+dd = dd) THEN  GOTO 2
         END;
         m := n;
2:         IF (m <> l) THEN BEGIN
            IF (iter = 30) THEN BEGIN
               writeln('pause in routine TQLI');
               writeln('too many iterations'); readln
            END;
            iter := iter+1;
            g := (d[l+1]-d[l])/(2.0*e[l]);
            r := sqrt(sqr(g)+1.0);
            g := d[m]-d[l]+e[l]/(g+sign(r,g));
            s := 1.0;
            c := 1.0;
            p := 0.0;
            FOR i := m-1 DOWNTO l DO BEGIN
               f := s*e[i];
               b := c*e[i];
               IF (abs(f) >= abs(g)) THEN BEGIN
                  c := g/f;
                  r := sqrt(sqr(c)+1.0);
                  e[i+1] := f*r;
                  s := 1.0/r;
                  c := c*s
               END ELSE BEGIN
                  s := f/g;
                  r := sqrt(sqr(s)+1.0);
                  e[i+1] := g*r;
                  c := 1.0/r;
                  s := s*c
               END;
               g := d[i+1]-p;
               r := (d[i]-g)*s+2.0*c*b;
               p := s*r;
               d[i+1] := g+p;
               g := c*r-b;
            (* Next loop can be omitted if eigenvectors not wanted *)
               FOR k := 1 TO n DO BEGIN
                  f := z[k,i+1];
                  z[k,i+1] := s*z[k,i]+c*f;
                  z[k,i] := c*z[k,i]-s*f
               END
            END;
            d[l] := d[l]-p;
            e[l] := g;
            e[m] := 0.0;
            GOTO 1
         END
      END
 END;


But I would like to know how to remove goto
and reduce space complexity
We need only one entry for each eigenvector to calculate weight

 

Laksen

  • Hero Member
  • *****
  • Posts: 785
    • J-Software
Re: Gaussian Quadratures
« Reply #1 on: July 12, 2024, 02:46:15 pm »
I cleaned it up a bit and removed the goto's, changed datatypes, raised exception instead of readln.

Not sure what you mean with space complexity? Do you need the eigenvectors or not?

Code: Pascal  [Select][+][-]
  1. uses
  2.   SysUtils;
  3.  
  4. type
  5.   glnp = array of double;
  6.   glnpnp = array of glnp;
  7.  
  8. function sign(a, b: double): double;
  9. begin
  10.   if (b < 0) then sign := -abs(a)
  11.   else
  12.     sign := abs(a);
  13. end;
  14.  
  15. // Programs using routine TQLI must define the types where np is the physical dimension of the matrix to be analyzed.
  16. procedure tqli(var d, e: glnp; n: longint; var z: glnpnp);
  17. var
  18.   i, m, l, iter, k: longint;
  19.   s, r, p, g, f, dd, c, b: double;
  20. begin
  21.   for i := 2 to n do
  22.     e[i - 1] := e[i];
  23.   e[n] := 0.0;
  24.  
  25.   for l := 1 to n do
  26.   begin
  27.     for iter := 0 to 30 do
  28.     begin
  29.       m := n;
  30.       for i := l to n - 1 do
  31.       begin
  32.         dd := abs(d[i]) + abs(d[i + 1]);
  33.         if (abs(e[i]) + dd = dd) then
  34.         begin
  35.           m := i;
  36.           break;
  37.         end;
  38.       end;
  39.       if m <> l then
  40.       begin
  41.         if (iter = 30) then
  42.           raise Exception.Create('TQLI: Too many iterations');
  43.         g := (d[l + 1] - d[l]) / (2.0 * e[l]);
  44.         r := sqrt(sqr(g) + 1.0);
  45.         g := d[m] - d[l] + e[l] / (g + sign(r, g));
  46.         s := 1.0;
  47.         c := 1.0;
  48.         p := 0.0;
  49.         for i := m - 1 downto l do
  50.         begin
  51.           f := s * e[i];
  52.           b := c * e[i];
  53.           if abs(f) >= abs(g) then
  54.           begin
  55.             c := g / f;
  56.             r := sqrt(sqr(c) + 1.0);
  57.             e[i + 1] := f * r;
  58.             s := 1.0 / r;
  59.             c := c * s;
  60.           end
  61.           else
  62.           begin
  63.             s := f / g;
  64.             r := sqrt(sqr(s) + 1.0);
  65.             e[i + 1] := g * r;
  66.             c := 1.0 / r;
  67.             s := s * c;
  68.           end;
  69.           g := d[i + 1] - p;
  70.           r := (d[i] - g) * s + 2.0 * c * b;
  71.           p := s * r;
  72.           d[i + 1] := g + p;
  73.           g := c * r - b;
  74.           // Next loop can be omitted if eigenvectors not wanted
  75.           for k := 1 to n do
  76.           begin
  77.             f := z[k, i + 1];
  78.             z[k, i + 1] := s * z[k, i] + c * f;
  79.             z[k, i] := c * z[k, i] - s * f;
  80.           end;
  81.         end;
  82.         d[l] := d[l] - p;
  83.         e[l] := g;
  84.         e[m] := 0.0;
  85.         continue;
  86.       end
  87.       else
  88.         break;
  89.     end;
  90.   end;
  91. end;
  92.  

Thaddy

  • Hero Member
  • *****
  • Posts: 16196
  • Censorship about opinions does not belong here.
Re: Gaussian Quadratures
« Reply #2 on: July 12, 2024, 05:44:02 pm »
Somehow I know that code  ;) but I wrote it in TB, not Pascal...Can be wrong, but somehow?
anyway, laksen, thanx for fixing this.
If I smell bad code it usually is bad code and that includes my own code.

user07022024

  • New Member
  • *
  • Posts: 13
Re: Gaussian Quadratures
« Reply #3 on: July 12, 2024, 05:51:02 pm »
According to Golub-Welsch eigenvalues of this type of matrix can be used as nodes for Gaussian Quadratures
and for calculating weights it is enough to calculate first entry of each eigenvector

Thaddy

  • Hero Member
  • *****
  • Posts: 16196
  • Censorship about opinions does not belong here.
Re: Gaussian Quadratures
« Reply #4 on: July 12, 2024, 06:39:33 pm »
That is correct. That is why laksen wrote his comment.
If I smell bad code it usually is bad code and that includes my own code.

Laksen

  • Hero Member
  • *****
  • Posts: 785
    • J-Software
Re: Gaussian Quadratures
« Reply #5 on: July 12, 2024, 07:25:30 pm »
It is only doing column operations on the z matrix. So if you just need the first component then it should be easy to reduce it to the first row of the matrix

VTwin

  • Hero Member
  • *****
  • Posts: 1224
  • Former Turbo Pascal 3 user
Re: Gaussian Quadratures
« Reply #6 on: July 12, 2024, 09:15:06 pm »
Just a note, I would be careful using Numerical Recipes code, even the old 1986 Pascal "shareware diskette", as it does not use a standard license.
“Talk is cheap. Show me the code.” -Linus Torvalds

Free Pascal Compiler 3.2.2
macOS 14.5: Lazarus 3.4 (64 bit Cocoa M1)
Ubuntu 18.04.3: Lazarus 3.4 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 3.4 (64 bit on VBox)

Laksen

  • Hero Member
  • *****
  • Posts: 785
    • J-Software
Re: Gaussian Quadratures
« Reply #7 on: July 12, 2024, 09:48:50 pm »
The paragraph in the start of the book is fairly clear (at least the 1992 reprint I have).
If you own the book you can make use of it for a program :)

The newer licenses seem a bit more cumbersome: https://numerical.recipes/licenseinfo.html

VTwin

  • Hero Member
  • *****
  • Posts: 1224
  • Former Turbo Pascal 3 user
Re: Gaussian Quadratures
« Reply #8 on: July 12, 2024, 10:53:24 pm »
The file "nrpas13.zip" (I can't find an active link now) includes readme docs about using the 1986 code. My copy of the Pascal book includes a long section on usage, "Legal Matters".

I used some of these routines for many years, but eventually removed them for open source alternatives with a clear license.

Just my two cents, use as you please.


 
“Talk is cheap. Show me the code.” -Linus Torvalds

Free Pascal Compiler 3.2.2
macOS 14.5: Lazarus 3.4 (64 bit Cocoa M1)
Ubuntu 18.04.3: Lazarus 3.4 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 3.4 (64 bit on VBox)

Curt Carpenter

  • Hero Member
  • *****
  • Posts: 566
Re: Gaussian Quadratures
« Reply #9 on: July 12, 2024, 10:56:28 pm »
There's always a first time I guess, but does anybody know if the authors or publisher of Numerical Recipes has ever sued anyone in the last 38 years?  (Or if anybody has sued them for bugs...) 
« Last Edit: July 12, 2024, 10:58:40 pm by Curt Carpenter »

VTwin

  • Hero Member
  • *****
  • Posts: 1224
  • Former Turbo Pascal 3 user
Re: Gaussian Quadratures
« Reply #10 on: July 12, 2024, 11:03:25 pm »
Ha. Probably not, I try to follow the open source rules though.
“Talk is cheap. Show me the code.” -Linus Torvalds

Free Pascal Compiler 3.2.2
macOS 14.5: Lazarus 3.4 (64 bit Cocoa M1)
Ubuntu 18.04.3: Lazarus 3.4 (64 bit on VBox)
Windows 7 Pro SP1: Lazarus 3.4 (64 bit on VBox)

Curt Carpenter

  • Hero Member
  • *****
  • Posts: 566
Re: Gaussian Quadratures
« Reply #11 on: July 13, 2024, 12:14:57 am »
And, sadly, your legal department probably owes you their thanks for being so conscientious.  (I live in Texas, the home of some of the world's most notorious "patent pirates" that make their living off things like license violations.)

jamie

  • Hero Member
  • *****
  • Posts: 6735
Re: Gaussian Quadratures
« Reply #12 on: July 13, 2024, 05:52:33 pm »
And, sadly, your legal department probably owes you their thanks for being so conscientious.  (I live in Texas, the home of some of the world's most notorious "patent pirates" that make their living off things like license violations.)

Ah, you mean ambulance chasers (Lawers). :o
The only true wisdom is knowing you know nothing

 

TinyPortal © 2005-2018