Forum > Beginners

(1/3) > >>

user07022024:

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

--- Code: ---
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');
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;

--- End code ---

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:
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  [+][-]window.onload = function(){var x1 = document.getElementById("main_content_section"); if (x1) { var x = document.getElementsByClassName("geshi");for (var i = 0; i < x.length; i++) { x[i].style.maxHeight='none'; x[i].style.height = Math.min(x[i].clientHeight+15,306)+'px'; x[i].style.resize = "vertical";}};} ---uses  SysUtils; type  glnp = array of double;  glnpnp = array of glnp; function sign(a, b: double): double;begin  if (b < 0) then sign := -abs(a)  else    sign := abs(a);end; // Programs using routine TQLI must define the types where np is the physical dimension of the matrix to be analyzed.procedure tqli(var d, e: glnp; n: longint; var z: glnpnp);var  i, m, l, iter, k: longint;  s, r, p, g, f, dd, c, b: double;begin  for i := 2 to n do    e[i - 1] := e[i];  e[n] := 0.0;   for l := 1 to n do  begin    for iter := 0 to 30 do    begin      m := n;      for i := l to n - 1 do      begin        dd := abs(d[i]) + abs(d[i + 1]);        if (abs(e[i]) + dd = dd) then        begin          m := i;          break;        end;      end;      if m <> l then      begin        if (iter = 30) then          raise Exception.Create('TQLI: Too many iterations');        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;        continue;      end      else        break;    end;  end;end;

Somehow I know that code  ;) but I wrote it in TB, not Pascal...Can be wrong, but somehow?
anyway, laksen, thanx for fixing this.

user07022024:
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