Сингулярное разложение в линейной задаче метода наименьших квадратов

Дипломная работа - Математика и статистика

Другие дипломы по предмету Математика и статистика

Скачать Бесплатно!
Для того чтобы скачать эту работу.
1. Пожалуйста введите слова с картинки:

2. И нажмите на эту кнопку.
закрыть



Лоусон Ч., Хенсон Р. Численное решение задач наименьших квадратов. М.: Статистика, 1979, 447с

  • Марчук Г.И. Методы вычислительной математики. М.: Наука, 1980
  • Мэйндоналд Дж. Вычислительные алгоритмы в прикладной статистике. М.: Финансы и статистика, 1988, 350с
  • Стренг Г. Линейная алгебра и ее применения. М.: Мир, 1980, 454с
  • Уилкинсон Дж., Райнш К. Справочник алгоритмов на языке АЛГОЛ. Линейная алгебра, М.: Машиностроение, 1976, 390с
  • Фаддеев Д.К., Фаддеева В.Н. Вычислительные методы линейной алгебры. -М.: Физматгиз, 1963, 536с.
  • Форсайт Дж., Малькольм М., Моулер К. Машинные методы математических вычислений. М.: Мир, 1980, 279с
  • Харебов К.С. Компьютерные методы решения задачи наименьших квадратов и проблемы собственных значений. Владикавказ.: Изд-во СОГУ, 1995, 76 с.
  • ПРИЛОЖЕНИЕ 1. Исходные тексты программы

    REAL A(3,3), U(3,3), V(3,3), SIGMA(3), WORK(3),Y(3),C(3),Y0(3)

    INTEGER I,IERR, J, M, N, NM

    OPEN (6,FILE="SVD.OUT",STATUS="UNKNOWN",FORM="FORMATTED")

    OPEN (5,FILE= "SVD.IN",STATUS="UNKNOWN",FORM="FORMATTED")

    140 FORMAT(3I5)

    150 FORMAT(4E15.7)

    READ(5,140) NM,M,N

    DO 131 I=1,M

    READ(5,150) (A(I,J),J=1,N)

    131 CONTINUE

    READ (5,150) (Y(I),I=1,M)

    CALL SVD(NM,M,N,A,SIGMA,.TRUE.,U,.TRUE.,V,IERR,WORK)

    IF(IERR.NE.0) WRITE (6,2) IERR

    2FORMAT(15H TROUBLE.IERR=,I4)

    WRITE(6,120)

    120 FORMAT(/МАТРИЦА А)

    DO 121 I=1,M

    WRITE(6,130) (A(I,J),J=1,N)

    130 FORMAT(8E15.7)

    121 CONTINUE

    WRITE (6,160) (Y(I),I=1,N)

    160 FORMAT(/ПРАВЫЕ ЧАСТИ/8E15.7)

    210 FORMAT(/СИНГУЛЯРНЫЕ ЧИСЛА)

    WRITE(6,210)

    DO 3 J=1,N

    WRITE(6,6) SIGMA(J)

    3CONTINUE

    SMA=SIGMA(1)

    SMI=SIGMA(1)

    DO 211 J=2,N

    IF(SIGMA(J).GT.SMA) SMA=SIGMA(J)

    IF(SIGMA(J).LT.SMI.AND.SIGMA(J).GT.0.) SMI=SIGMA(J)

    211 CONTINUE

    OBU=SMA/SMI

    230 FORMAT(/ЧИСЛО ОБУСЛОВЛЕННОСТИ=,E15.7)

    WRITE(6,230) OBU

    SIGMA1=0.

    DO 30 J=1,N

    IF(SIGMA(J) .GT. SIGMA1) SIGMA1=SIGMA(J)

    C(J)=0.

    30 CONTINUE

    TAU=SIGMA1*0.1E-6

    DO 60 J=1,N

    IF(SIGMA(J).LE.TAU) GO TO 60

    S=0.

    DO 40 I=1,N

    S=S+U(I,J)*Y(I)

    40 CONTINUE

    S=S/SIGMA(J)

    DO 50 I=1,N

    C(I)=C(I) + S*V(I,J)

    50 CONTINUE

    60 CONTINUE

    write (6,560)

    WRITE (6,6) (C(I),I=1,3)

    DO 322 J=1,N

    SS=0.

    DO 321 I=1,M

    321 SS=A(J,I)*C(I)+SS

    322 Y0(J)=SS

    write (6,570)

    WRITE (6,6) (Y0(I),I=1,3)

    CWRITE(6,7)

    CDO 4 I=1,M

    CWRITE(6,6) (U(I,J),J=1,N)

    C4CONTINUE

    CWRITE(6,7)

    CDO 5 I=1,N

    CWRITE(6,6) (V(I,J),J=1,N)

    C5CONTINUE

    6FORMAT(3E15.7)

    560format(2x,roots)

    570format(2x,right)

    7FORMAT(1H )

    STOP

    E N D

    SUBROUTINE SVD(NM,M,N,A,W,MATU,U,MATV,V,IERR,RV1)

    REAL A(NM,N),W(N),U(NM,N),V(NM,N),RV1(N)

    LOGICAL MATU,MATV

    IERR=0

    DO 100 I=1,M

    DO 100J=1,N

    U(I,J)=A(I,J)

    100 CONTINUE

    G=0.0

    SCALE=0.0

    ANORM=0.0

    DO 300 I=1,N

    L=I+1

    RV1(I)=SCALE*G

    G=0.0

    S=0.0

    SCALE=0.0

    IF(I.GT.M) GO TO 210

    DO 120 K=I,M

    120SCALE=SCALE+ABS(U(K,I))

    IF(SCALE.EQ.0.0) GO TO 210

    DO 130 K=I,M

    U(K,I)=U(K,I)/SCALE

    S=S+U(K,I)**2

    130 CONTINUE

    F=U(I,I)

    G=-SIGN(SQRT(S),F)

    H=F*G-S

    U(I,I)=F-G

    IF(I.EQ.N) GO TO 190

    DO 150 J=L,N

    S=0.0

    DO 140 K=I,M

    140S=S+U(K,I)*U(K,J)

    F=S/H

    DO 150 K=I,M

    U(K,J)=U(K,J)+F*U(K,I)

    150CONTINUE

    190 DO 200 K=I,M

    200U(K,I)=SCALE*U(K,I)

    210W(I)=SCALE*G

    G=0.0

    S=0.0

    SCALE=0.0

    IF(I.GT.M.OR.I.EQ.N) GO TO 290

    DO 220 K=L,N

    220SCALE=SCALE+ABS(U(I,K))

    IF(SCALE.EQ.0.0) GO TO 290

    DO 230 K=L,N

    U(I,K)=U(I,K)/SCALE

    S=S+U(I,K)**2

    230CONTINUE

    F=U(I,L)

    G=-SIGN(SQRT(S),F)

    H=F*G-S

    U(I,L)=F-G

    DO 240 K=L,N

    240RV1(K)=U(I,K)/H

    IF(I.EQ.M) GO TO 270

    DO 260 J=L,M

    S=0.0

    DO 250 K=L,N

    250S=S+U(J,K)*U(I,K)

    DO 260 K=L,N

    U(J,K)=U(J,K)+S*RV1(K)

    260CONTINUE

    270DO 280 K=L,N

    280U(I,K)=SCALE*U(I,K)

    290ANORM=AMAX1(ANORM,ABS(W(I))+ABS(RV1(I)))

    300CONTINUE

    IF(.NOT.MATV) GO TO 410

    DO 400 II=1,N

    I=N+1-II

    IF(I.EQ.N) GO TO 390

    IF(G.EQ.0.0) GO TO 360

    DO 320 J=L,N

    320V(J,I)=(U(I,J)/U(I,L))/G

    DO 350 J=L,N

    S=0.0

    DO 340 K=L,N

    340S=S+U(I,K)*V(K,J)

    DO 350 K=L,N

    V(K,J)=V(K,J)+S*V(K,I)

    350CONTINUE

    360DO 380 J=L,N

    V(I,J)=0.0

    V(J,I)=0.0

    380CONTINUE

    390V(I,I)=1.0

    G=RV1(I)

    L=I

    400CONTINUE

    410IF(.NOT.MATU) GO TO 510

    MN=N

    IF(M.LT.N) MN=M

    DO 500 II=1,MN

    I=MN+1-II

    L=I+1

    G=W(I)

    IF(I.EQ.N) GO TO 430

    DO 420 J=L,N

    420U(I,J)=0.0

    430IF(G.EQ.0.0) GO TO 475

    IF(I.EQ.MN) GO TO 460

    DO 450 J=L,N

    S=0.0

    DO 440 K=L,M

    440S=S+U(K,I)*U(K,J)

    F=(S/U(I,I))/G

    DO 450 K=I,M

    U(K,J)=U(K,J)+F*U(K,I)

    450CONTINUE

    460DO 470 J=I,M

    470U(J,I)=U(J,I)/G

    GO TO 490

    475DO 480 J=I,M

    480U(J,I)=0.0

    490U(I,I)=U(I,I)+1.0

    500CONTINUE

    510DO 700 KK=1,N

    K1=N-KK

    K=K1+1

    ITS=0

    520DO 530 LL=1,K

    L1=K-LL

    L=L1+1

    IF(ABS(RV1(L))+ANORM.EQ.ANORM) GO TO 565

    IF(ABS(W(L1))+ANORM.EQ.ANORM) GO TO 540

    530CONTINUE

    540C=0.0

    S=1.0

    DO 560 I=L,K

    F=S*RV1(I)

    RV1(I)=C*RV1(I)

    IF(ABS(F)+ANORM.EQ.ANORM) GO TO 565

    G=W(I)

    H=SQRT(F*F+G*G)

    W(I)=H

    C=G/H

    S=-F/H

    IF(.NOT.MATU) GO TO 560

    DO 550 J=1,M

    Y=U(J,L1)

    Z=U(J,I)

    U(J,L1)=Y*C+Z*S

    U(J,I)=-Y*S+Z*C

    550CONTINUE

    560CONTINUE

    565Z=W(K)

    IF(L.EQ.K) GO TO 650

    IF(ITS.EQ.30) GO TO 1000

    ITS=ITS+1

    X=W(L)

    Y=W(K1)

    G=RV1(K1)

    H=RV1(K)

    F=((Y-Z)*(Y+Z)+(G-H)*(G+H))/(2.0*H*Y)

    G=SQRT(F*F+1.0)

    F=((X-Z)*(X+Z)+H*(Y/(F+SIGN(G,F))-H))/X

    C=1.0

    S=1.0

    DO 600 I1=L,K1

    I=I1+1

    G=RV1(I)

    Y=W(I)

    H=S*G

    G=C*G

    Z=SQRT(F*F+H*H)

    RV1(I1)=Z

    C=F/Z

    S=H/Z

    F=X*C+G*S

    G=-X*S+G*C

    H=Y*S

    Y=Y*C

    IF(.NOT.MATV) GO TO 575

    DO 570 J=1,N

    X=V(J,I1)

    Z=V(J,I)

    V(J,I1)=X*C+Z*S

    V(J,I)=-X*S+Z*C

    570CONTINUE

    575Z=SQRT(F*F+H*H)

    W(I1)=Z

    IF(Z.EQ.0.0) GO TO 580

    C=F/Z

    S=H/Z

    580F=C*G+S*Y

    X=-S*G+C*Y

    IF(.NOT.MATU) GO TO 600

    DO 590 J=1,M

    Y=U(J,I1)

    Z=U(J,I)

    U(J,I1)=Y*C+Z*S

    U(J,I)=-Y*S+Z*C

    590CONTINUE

    600CONTINUE

    RV1(L)=0.0

    RV1(K)=F

    W(K)=X

    GO TO 520

    650IF(Z.GE.0.0) GO TO 700

    W(K)=-Z

    IF(.NOT.MATV) GO TO 700

    DO 690 J=1,N

    690V(J,K)=-V(J,K)

    700CONTINUE

    GO TO 1001

    1000IERR=K

    1001RETURN

    E N D

    ПРИЛОЖЕНИЕ 2. контрольный пример

    Входные данные

    (матрица изначально сингулярна первая строка равна сумме второй и третьей с обратным знаком)

    3 3 3

    .3200000E 02 .1400000E 02 .7400000E 02

    -0.2400000E 02 -0.1000000E 02 -0.5700000E 02

    -0.8000000E 01 -0.4000000E 01 -0.1700000E 02

    -0.1400000E 02 0.1300000E 02 0.1000000E 01

     

    Полученный результат

    МАТРИЦА А

    .3200000E+02 .1400000E+02 .7400000E+02

    -.2400000E+02 -.1000000E+02 -.5700000E+02

    -.8000000E+01 -.4000000E+01 -.1700000E+02

     

    ПРАВЫЕ ЧАСТИ

    -.1400000E+02 .1300000E+02

    s