|
Главная -> Справочник по алгоритмам 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 [38] 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 Корни полиномов Р (%) могут быть вычислены по программам, описанным в § 4.5. Максимальное собственное значение матрицы (действительное) и принадлежащий ему собственный вектор могут вычисляться степенным методом по следующему алгоритму. 1. По формуле Y>=Ay-\ где ац ai2 021 022 а,„ а2п Loni ап2 . вычисляется вектор Y (перед началом вычислений задается произвольный исходный вектор Y"). 2. Находятся приближения к наибольшему собственному значению 3. Вектор Y- нормируется, т. е. вместо него получается вектор 4. Проверяется выполнение условия 1р.с-Hr-ll <Е, где Е - заданная погрешность вычисления максимального по модулю собственного значения. Если это условие выполняется, считаем lc максимальным по модулю собственным значением, а вектор f* - принадлежащим ему собственным вектором. Если данное условие не выполняется, идем к п. 1, приняв за У<~" вектор УЧ Программа 4.58. твердого материала со всех сторон действуют силы, которые в матричной форме записываются в виде (все элементы матрицы нормированы, реальные усилия получаем их умножением на 10" Н/м) 10 5 6 5 20 4 6 4 30 При N = 3 и £ = £=1-10"" находим Цыакс = = /-«акс = 33,70906694. Кроме того, вычисляется принадлежащий ему вектор: yi = = 0,3408143922, j/2 = 0,4160822399 и уз=1. Все действительные собственные значения симметричной матрицы {илк заданное их число) могут определяться методом скалярных произведений, реализованным с помощью следующего алгоритма. 1. Вычисляем вектор К**=Л,У<*~где k - номер итерации. Л, - матрица А. Находим начальные приближения к максимальному собственному значению fi*=(i/*". {/*")/({/-". i/*")- 3. Нормируем вектор У**, т. е. заменяем его вектором 4. Проверяем выполнение неравенства Если оно выполняется, считаем ц* наибольшим собственным значением матрицы Ai, а Р* - принадлежащим ему собственным вектором. Если данное условие не выполняется, повторяем вычисления с п. 1, приняв за вектор f. 5. Если задано число искомых собственных значений г=1, 2, k, то для каждого i повторяем вычисления по п. 1-4, преобразовав матрицу Ai по формуле Ai=Ai-.,-%i [Ui, и/}, i=l, 2.....k-\. где [Ui, Ui] - произведение вектора-столбца на вектор-строку. 1в PRINTВЫЧИСЛЕНИЕ МАКСИМАЛЬНОГО СОБСТВЕННОГО ЗНАЧЕНИЯ 15 PRINTMATPHUH И ПРИНАДЛЕМАИЕГО ЕМУ НОРМИРОВАННОГО ВЕКТОРА 20 ШРиТВВЕДИТЕ И,Е N/E:DIM A<M/N):FOR 1=1 ТО H:FOR J=l ТО N 30 PRINTBBEAHTE А12.0!I/J:INPUT A<I,J)iNEXT JiNEXT I 40 FOR I»l TO N: PRINT 12.01ВВЕДИТЕ УвСГ) 50 INPUT A<e,l)sNEXT IiiLETU=e 66 LETL=e:LETC»=0«LETR=0:FOR J=l TO HsLETS=0 70 FOR 1-1 TO N:LETS=S+A<bJ)a«A<e/I):NEXT I 80 LETA<J/e)=S:IF ABS<S)>R THEN LETR=ABS<S) 90 IF S«8 THEN 120 lee IF A<0/J)=0 THEN 120 110 LETC=C+1«LETL*L+A<J*0)/A<0/J) 120 NEXT J:LETL=L/t; 130 FOR J=l TO NsLETA<e/J)=A<J/0)/R:NEXT J 148 IF ABS<L-U)>«E THEN LETU=L!60T0 66 150 PRINT«F1.9«МАКСИМАЛЬНОЕ COBCTBEHHOE ЗНАЧЕНИЕ L=L 168 PRINTBEKTOPiFOR I-l TO N 170 PRlNT!2.e!VI»!F1.9!A<0/l)SHEXT IsEND Пример. Для контроля этой программы воспользуемся примером из [41J: на куб Этот метод применим, если все собственные значения матрицы А действительны. а наибольшее собственное значение не кратно. Перед началом вычислений вводится матрица А, задается погрешность вычислений £ = число искомых собственных значений и векторов К (1</СЛ/ = п) и начальный вектор К=и. Для экономии памяти ПЭВМ симметричная матрица задается своим нижним треугольником построчно и занимает одномерный массив с объемом элементов, равным N(N-{-l)/2. Программа 4.59. 2. Вычисляем первую преграду ,7=1 3. Находим внедиагональный элемент а,,, по модулю превосходящий а/, {k = 0, 1,2,...). Если его нет, выполняем п. 5 алгоритма, иначе идем к п. 4. 10 PRINTВЫЧИСЛЕНИЕ СОБСТВЕННЫХ ЗНАЧЕНИЙ И ПРИНАДЛЕЖАЩИХ Ж 15 PRINTBEKTOPOB ДЛЯ СИММЕТРИЧНОЙ МАТРИЦЫ МЕТОДОМ СКАЛЯРНЫХ 20 PRINT ПРОИЗВЕДЕНИЙ С ИСЧЕРПЫВАНИЕМ 25 1НРиТВВЕДИТЕ N/E/K N»E/K!LETU=H*<N+l)/2;LETU=W+N:mm A<U+N) 30 LETH=0:FOR 1=1 ТО N:FOR J=l TO l!LETH=H+l 40 PRINT12.0!ВВЕДИТЕ ftI*J! INPUT A<H>; 50 NEXT JsNEXT UFOR B=l TO К 60 FOR 1=1 TO NSLETA<U+I)=1!NEXT i:LETF=« 70 FOR 1=1 TO NSLETA<U+IJ=0:FOR J=l TO I 80 LETA<U+I)=A<U+I)+aa*a-l>2+J)iiiA<U+J):NEXT J 90 NEXT irFOR J=l TO N-lsFOR I=J+1 TO N 100 LETA<U+J)=A<U+J>+A(l*<I-l)+J)*ft(U+I) 110 NEXT ISNEXT JsLETP=0sLETR=0 120 FOR 1=1 TO NsLETP=P+A<U+l)*ft<U+I) 130 LETR=R+A<W+I>i*A<U+l):NEXT I 140 LETL=P/RsLETC=SQR<P)sLETP=0 150 FOR 1=1 TO NsLETA<U+I>=ft<U+I) 160 IF ABS«:A<U+1)>>=P THEN LETP=ABS<A«:W+I)) 170 NEXT IS IF ABS<P-F)>E THEN LETF=Ps60T0 70 100 PRINT СОБСТВЕННОЕ ЗНАЧЕНИЕ L И ЕГО ВЕКТОР V 190 PRINT!2.0!LI=!F1.9!L 200 FOR 1=1 TO HsPRINT!2.0!VI=!F1.9!A<:U+n 210 NEXT niF B=K THEH 248 220 FOR 1=1 TO NSFOR J=l TO IsLETB=H*<I-l)/2+J 230 LETA«:D)=A<D)-L«A(U+I)5I!A<U+J)sNEXT JSNEXT I 240 NEXT BSEND Пример. Для матрицы из примера к программе 4.58 получаем: 1) /. = 33,70917846. {/1=0,300158658, {/2 = 0,3664623805, Уз = 0,8806872904; 2) /.= 19,14906125, {/, =0,1967307594, {/2 = 0,8796183084, j3= 0,4330919532 и 3) /. = 7,141760291, у, =0,933376197, {/2=-0,3032741433, Уз =-0,19192099959. Для отыскания всех собственных значений матрицы применяется также метод Якоби с преградами. Его суть заключается в проведении цепочки преобразований подобия, в ходе которых получается некоторая диагональная матрица Л**: Л=7Л7, имеющая те же собственные значения, что и матрица Л. Матрица Т является произведением всех матриц 7",,, где Tq - элементарные матрицы вращения, имеющие вид 4. .Анализируем элемент с,-,-. Для этого вычисляем если у = 0. Тг.= ...S -s:-.c c4s==i. Алгоритм реализации метода Якоби с преградами следующий. 1. Задаем Тт = Е. где £ -единичная матрица порядка п. У=(а„ -а,-,)/2, ;,=/ , I -sign (у) aiiafj+y, если уф О, 5=х/д/2 (1-Vl-W. C = Vl-S=. Затем преобразуем г-й и / й столбцы матрицы Л* по формулам B, = CAi-SAi, Bi = SAi + CAi и заменяем (-й и /-й столбцы матрицы Л столбцами Bi и В,-. В результате на месте матрицы Л получаем матрицу В. Далее (-Ю и ;-ю строки матрицы В преобразуем по формулам C, = СЛ,-5Л,-, Ci = SA, + CAi и заменяем t-ю и . ;-ю строки матрицы В строками С; и Cj. При этом вместо матрицы В получаем матрицу Л**. Столбцы i и / матрицы Тц заменяем на столбцы В, и В,. 5. Находим новую преграду «*+,=«*/«(="+"+> и повторяем вычисления с п. 1, до тех пор, пока все внедиагональные элементы не станут меньше числа Р = еяс, где е - заданная погрешность вычислений. В результате собст- венные значения оказываются диагональ- Пример. Для матрицы А из примера ными элементами матрицы А. к программе 4.58 вычисления по программе Ввиду симметрии матрицы она задается 4.60 дают Z,i =7,14176. Z,2= 19,14906 и построчно нижним треугольником и зани- Z,3 = 33,709I7 (/V = 3, г = £= I 10~). мает одномерный массив с числом элементов Программы 4.57-4.60 получены перево- N(N+\)/2. . дом программ для ЭВМ «Мир-2» с языка Программа 4.60. аналитик [24] на язык бейсик. 10 PRINTВЫЧИСЛЕНИЕ ВСЕХ СОБСТВЕЖЫХ ЗНАЧЕНИЙ МАТРИЦЫ 15 PRINT МЕТОДОМ ЯКОБИ С ПРЕГРАДАМИ 20 INPUTBBEAHTE ПОРЯДОК МАТРИЦЫ N=N:DIM A<N*(N+1)2) 25 INPUTBBEAHTE ПОГРЕШНОСТЬ ВЫЧИСЛЕНИЯ Е=Е 30 LETH=0:FOR 1=1 ТО NsFOR J=l ТО I 40 LETH=H+1!PRINT!3.0!ВВЕДИТЕ АI/JsINPUT А(Н) 50 NEXT JSNEXT I:LETH=0:FOR P=2 TO N 60 FOR Q=l TO P-1:LETH=H+2*A<:P*<:P-1)2+S>2 70 NEXT QSNEXT p!LETR=SeR<:H>!LETA=Ei«R/H 80 LETR=RN Э0 LETB=e!FOR 0=2 TO NsFOR P=l TO G-1 188 IF ABStAte*tG-l>2+P)><R THEN 275 110 LETK=P*<P-1 > 2+Ps LETL=ei«( Q-l >2+P 120 LETM=S*(e-l)2+QsLETB=l!LETU=AtK>sLETW=AtL) 130 LETF=A(M>sLETV=tU-F).2sLETZ=ViIF V=8 THEN LETX=-l!60TO 158 140 LETX=-S6NtZ>3tiWSQRtl«+V*V> 158 LETS=X/SQRt2*<l+SQRa-X*X>>)sLETC=SeR<l-S*S) 160 FOR 1=1 TO NsIF I<=P THEN LETV=P*tP-l)2+Is60T0 188 178 -LETV=I*a-U/2+P 138 IF I<=G THEN LETX=e*(Q-l>/2+ls60T0 288 1Э0 LETX=I*a-l >2+U 200 IF I=G THEN 220 210 LETIi=AtV>iKC-A<:x>*S 228 LETAtX>=A(V>*S+AtX>*CsIF 1=0 THEN 248 238 L£TA<V)=D 240 NEXT I 258 LETAtK>=U*C*C+F*S*S-2!iiW*C*S 260 LETAtM)=U*S*S+F*C*C+2*W*C*S 278 LETAt L >=tU-F)жЗжС+Ы* t C*C-S*S) 275 NEXT PsNEXT Q 288 IF B=l THEN 98 290 IF R>A THEN LETR=RN!60T0 80 308 FOR 1=1 TO NSLETH=Aa»a-l).2+I) 318 PRINT!2.8!СОБСТВЕННОЕ ЗНАЧЕНИЕ LI=!F1.9!Н 328 NEXT ISEND 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 [38] 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 0.0164 |
|