Реализация QL-алгоритма в scilab

@xzkto · Регистрация: 20.06.2014

Author24 — интернет-сервис помощи студентам

Пытаюсь реализовать QL-алгоритм в scilab
Использовал текст программы на языке C из статьи http://www.cyberguru.ru/progra... age10.html

Matlab M

function[d,z]=ql(d,e,z)
    MAXITER=30;
    n=length(d);
    for l=1:1:n,
        iter=0, m=0,
        while 1==1 do
            for m=l:1:n-1,
                dd=abs(d(m))+abs(d(m+1)),
                if abs(e(m)+dd)==dd then break; end,
            end,
            if m<>l then
                iter=iter+1,
                if (iter)>=MAXITER break; end
                g=(d(l+1)-d(l))/2*e(l),
                r=sqrt(1+g^2),
                if g>=0 then g=g+abs(r),
                else g=g-abs(r), end
                g=d(m)-d(l)+e(l)/g,
                s=1,c=1,p=0,
                for i=m-1:-1:l,
                    f=s*e(i), b=c*e(i),
                    r=sqrt(f^2+g^2),
                    e(i+1)=r,
                    if r==0 then 
                        d(i+1)=d(i+1)-p, 
                        e(m)=0,
                        break;
                    end,
                    s=f/r, c=g/r, g=d(i+1)-p, r=(d(i)-g)*s+2*c*b, p=s*r, d(i+1)=g+r; g=c*r-b,
                    for k=1:1:n,
                        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,
                    if r==0 and i>=l then continue; end,
                    d(l)=d(l)-p, e(l)=g, e(m)=0,
                end,
            end,
            if l==m then break; end,
        end,
    end,
endfunction

Результат вызова:

Matlab M

d=[3.31,23.98,8.29,6.74];e=[5.34,9.05,0.98,0];z=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];
[d,z]=ql(d,e,z)
 !--error 21 
Неправильный индекс.
at line      14 of function ql called by :  
[d,z]=ql(d,e,z)

В программе на c# данная ошибка не возникала. Подскажите как исправить. Спасибо.

@Konimar · 20.06.2014, 19:24

Может так:

Matlab M

function[d,z]=ql(d,e,z)
    MAXITER=30;
    n=length(d);
    d=[d,0]    //!!!!!!!!!!!!!!!!
    for l=1:1:n,
        iter=0, 
        while 1==1 do
            for m=l:1:n-1,
                dd=abs(d(m))+abs(d(m+1)),
                if abs(e(m)+dd)==dd then break; end,
            end,
            if m<>l then
                iter=iter+1,
                if (iter)>=MAXITER break; end
                g=(d(l+1)-d(l))/2*e(l),
                r=sqrt(1+g^2),
                if g>=0 then g=g+abs(r),
                else g=g-abs(r), end
                g=d(m)-d(l)+e(l)/g,
                s=1,c=1,p=0,
                for i=m-1:-1:l,
                    f=s*e(i), b=c*e(i),
                    r=sqrt(f^2+g^2),
                    e(i+1)=r,
                    if r==0 then 
                        d(i+1)=d(i+1)-p, 
                        e(m)=0,
                        break;
                    end,
                    s=f/r, c=g/r, g=d(i+1)-p, r=(d(i)-g)*s+2*c*b, p=s*r, d(i+1)=g+r; g=c*r-b,
                    for k=1:1:n,
                        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,
                    if r==0 and i>=l then continue; end,
                    d(l)=d(l)-p, e(l)=g, e(m)=0,
                end,
            end,
            if l==m then break; end,
        end,
    end,
    d=d(1:n)       //!!!!!!!!! 
endfunction

Результат надо сравнить с результатом С

@xzkto · 21.06.2014, 02:15 **[ТС]**

Таким образом ошибка не возникает, но результаты совершенно другие.
Результаты вызова на C:

Matlab M

0,8883506
5,2234442
7,1653046
29,0429006
0,8068110 0,5237021 0,1997651 0,1867781
-0,3658826 0,1876545 0,1442238 0,9000640
0,4575097 -0,6979360 -0,3858370 0,3933193
-0,0766211 0,4510070 -0,8890575 0,0172826

Scilab:

Matlab M

[d,z]=ql(d,e,z)
 z  = 
    0.5587725    0.8083585    0.1852829    0.  
  - 0.3519487    0.0288402    0.9355749    0.  
    0.7509363  - 0.5879836    0.3006158    0.  
    0.           0.           0.           1.  
 d  =
 
  - 10.368253  - 9.7618808    42.927946    6.74

Результат при использовании встроеной функции:

Matlab M

@Konimar · 21.06.2014, 11:53

Правильного результата у меня не получается, но по сравнению с описанием на c нашел две ошибки :

Matlab M

1	g=(d(l+1)-d(l))/(2*e(l)) // - делить на 2e...

Matlab M

1	r=(d(i)-g)s+2cb, p=sr, d(i+1)=g+p; g=c*r-b, // было d(i+1)=g+r

@xzkto · 21.06.2014, 12:06 **[ТС]**

Вторую ошибку я нашел, а вот по поводу первой не понял, в чем проблема...
Немного переделал код в соответствии с http://beige.ucs.indiana.edu/B673/node38.html

Matlab M

function[d,z]=ql(d,e,z)
    MAXITER=30;
    n=length(d);
    d=[d,0],
    for l=1:1:n,
        iter=0,
        while 1==1 do
            for m=l:1:n-1,
                dd=abs(d(m))+abs(d(m+1)),
                if abs(e(m)+dd)==dd then break; end,
            end,
            if m==l then break; end,
            if (iter)>=MAXITER break; end
            iter=iter+1,
            g=(d(l+1)-d(l))/2*e(l),
            r=sqrt(1+g^2),
            if g>=0 then g=g+abs(r),
            else g=g-abs(r), end
            g=d(m)-d(l)+e(l)/g,
            s=1,c=1,p=0,
            for i=m-1:-1:l,
                f=s*e(i), b=c*e(i),
                r=sqrt(f^2+g^2),
                e(i+1)=r,
                if r==0 then 
                    d(i+1)=d(i+1)-p, 
                    e(m)=0,
                    break;
                end,
                s=f/r, c=g/r, g=d(i+1)-p, r=(d(i)-g)*s+2*c*b, p=s*r, d(i+1)=g+p; g=c*r-b,
                for k=1:1:n,
                    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,
            if l==m then break; end,
        end,
    end,
    d=d(1:n),
endfunction
 
d=[3.31,23.98,8.29,6.74];e=[5.34,9.05,0.98,0];
z=[1,0,0,0;0,1,0,0;0,0,1,0;0,0,0,1];
[d,z]=ql(d,e,z)
 z  = 
    0.8160482    0.5469301    0.1869033    0.  
  - 0.3649127    0.2367617    0.9004347    0.  
    0.4482233  - 0.8030015    0.3927907    0.  
    0.           0.           0.           1.  
 d  = 
    0.9221091    5.621644    29.036247    6.74

если говорить о собственных векторах, то результаты похожи на настоящие, правда точность очень плохая

@Konimar · 21.06.2014, 12:20

Скобок у Вас в строке 15 не хватает: (2*e(l))

Добавлено через 11 минут
Будет появляться деление на 0 -
против этого предлагаю переделать фрагмент:

Matlab M

if  m-1>=l then
                  g=(d(l+1)-d(l))/(2*e(l)),
                  r=sqrt(1+g^2),
                  if g>=0 then 
                      g=g+abs(r),
                  else 
                      g=g-abs(r), 
                  end
                  g=d(m)-d(l)+e(l)/g,
                  s=1,c=1,p=0,
                end;

Поскольку в следующем цикле только при l...m-1 значения обрабатываются

@xzkto · 21.06.2014, 12:21 **[ТС]**

Спасибо. Подставил скобки, в результате деление на нуль, т.к. вектор e=[5.34,9.05,0.98,0];
Снова вопрос, почему такой ошибки нету в C, при тех же входных данных?

@Konimar · 21.06.2014, 12:23

Я честно говоря с си знаком очень поверхностно, но там деление на 0 для действительных чисел не является ошибкой. А в последующем цикле эти данные (после деления н 0) не используются (по моему)

@xzkto · 21.06.2014, 18:33 **[ТС]**

Итак, дополнительное условие помогло избавиться от ошибки, собственные значения довольно близки к полученным встроенной функцией (но на С точность лучше), но собственные вектора считает не верно.
[d,z]=ql(d,e,z)

Matlab M

 z  =
     0.8160482    0.5469301    0.1869033    0.  
  - 0.3649127    0.2367617    0.9004347    0.  
    0.4482233  - 0.8030015    0.3927907    0.  
    0.           0.           0.           1.  
 d  = 
    0.9221091    5.621644    29.036247    6.7388739

@xzkto 0 / 0 / 0 Регистрация: 20.06.2014 Сообщений: 5
	21.06.2014, 12:21 [ТС]	7
	Спасибо. Подставил скобки, в результате деление на нуль, т.к. вектор e=[5.34,9.05,0.98,0]; Снова вопрос, почему такой ошибки нету в C, при тех же входных данных? 0

@Konimar 50 / 16 / 8 Регистрация: 01.10.2010 Сообщений: 77
	21.06.2014, 12:23	8
	Я честно говоря с си знаком очень поверхностно, но там деление на 0 для действительных чисел не является ошибкой. А в последующем цикле эти данные (после деления н 0) не используются (по моему) 0