Ai(n+i) сразу подставляются в последующие уравнения. Обычно (но не всегда) метод Зейделя дает лучщую сходимость.

Необходимость преобразования (4.16) и трудности в обеспечении быстрой сходимости ограничивают применение этих методов. Однако программная реализация их довольно проста (см. Приложение 5).

Метод Ньютона (или Ньютона - Раф-сона) является наиболее распространенным методом рещения системы уравнений (4.16). Он реализуется следующим алгоритмом.

1. Задаем абсолютную или относительную погрещиость в=£, число уравнений Л/, максимальное число итераций М и вектор начальных приближений Хю (с компонентами

XtO, Х20.....-Кда)-

2. Используя разложение Fi{Xi) в ряд Тейлора, формируем матрицу Якоби [dFi/dXi], необходимую для расчета приращений FXi) при малом изменении переменных. Матрица Якоби в развернутом виде записывается следующим образом:

дР,/дх, aF,/dxi ... dFi/dXf, dFi/дх, dF-i/dXi ... dF/dx

dFfi/dXt dF/dxi

дРц/дХц

Поскольку аналитическое дифференцирование Fi(Xi) в общем случае нежелательно, заменяем частные производные в матрице Якоби их приближенными конечно-разностными значениями

dFi FiXi + H)-FiXi)

3. Составляем и решаем систему линейных уравнений для малых приращений Хг.

где Hi - малое приращение Xi, nanpHiep Я. = е1Х,1.

Решение этой системы дает Дх), Ах2, ... т. е. ДЛ,.

4. Вычисляем уточненные значения

Х1(п+1) = Х(п) + ДХ1, Х2(„+ 1) = + ДХ2,

(п4-1) ~% (и) или в общем виде

Х- (л + 1) = •П +

5. Для всех AXi проверяем одно из условий: ДЛ,>е, ДХ,-/А:,->е.

Если оно выполняется, идем к п. 2, т. е. выполняем новую итерацию. Иначе считаем вектор Х(п-(-1) найденным решением.

Отметим, что решение системы нелинейных уравнений (4.16) можно представить и в виде Х,(„+ ,)=Хм- W(7)f,(A,-t„,), где Win) - обращенная матрица Якоби. Обращение матрицы Якоби осуществляется в ходе решения системы линейных уравнений для приращений ДХ, методом Гаусса. Программа 4.20.

16 PRINTРЕШЕНИЕ С1СТЕМЫ НЕЛИНЕЙНЫХ УРАВНЕНИЙ 26 PRINT . МОДИФИЦИРОВАННЫМ МЕТОДОМ НЫвТОНА

30 ШРиТЗАДАйТЕ ЧИСЛО УРАВНЕНИЙ N=NsriIM A<NN),B<N),X<N),F<N)

40 INPUTЗАДАЙТЕ МАКСИМАЛЬНОЕ ЧИСЛО ИТЕРАЦИЙ М=М

50 ШРИТЗАДАЙТЕ ОТНОСИТЕ ЛЬНУ» ПОГРЕШНОСТЬ E=E!LETS=e

66 FOR I-l ТО N!PRINT 12.0!ВВЕДИТЕ ХЧ<е)

76 IHPUT X<I>sHEXT I

80 60SUB 260:FOR 1=1 TO N!LETB<I)=-F<D!NEXT I

96 FOR J-1 TO N:LETX=X<d)sLETH=EiKABS<X>

166 LETX<J)=X+Hi60SUB 2601FOR 1=1 TO N

116 LETA<bJ)=<F<I)+B<I))/H!NEXT I:LETX<J>=X!NEXT J

120 LETS-S+ISIF S=M+1 THEN PRINTЧИСЛО ИТЕРАЦИЙ S=M:STOP

136 FOR 1=1 TO N-i:FOR J=I+1 TO N

140 leta<J,I>=-A<J,I>A<bI>:for K=I+1 TO N

150 L£TA<JK>=A<J,K>+A<J,niKA<IK):NEXT К

160 LETB<J)=B<J)+A<v),I)*B<I>!NEXT J!NEXT I

170 I.ETF<N>=B<N>/A<:N*N>!F0R I=N-1 to 1 STEP -1

186 LETH=B<I>»for J=l+1 TO N!LETH=H-F<J>iKA<bJ):NEXT J

190 LETF<I)=H-A<bI):NEXT I-.LETR=6

206 FOR 1=1 TO NsLETX<I>=X<I>+F<I)

210 IF ABS<F<I)X<I))>E THEN LETR=1

226 NEXT I:IF R=l THEN 86

230 PRINTРЕШЕНИЕ СИСТЕМЫ*

240 FOR 1=1 TO N:PRINT!2.6!XI!F1.9! = X(I)!NEXT I-245 PRINT!2.0!ЧИСЛО ИТЕРАЦИЙ S=S!STOP

250 REMПОДПРОГРАММА ВЫЧИСЛЕНИЯ F(I>=F<X<1>,X<2),...*Х<Н)> 266 LETF<l>=X<i)+3*L6T<X<l>)-X<2>*X<2) 270 ьеТР<2)=г»Х<1>»Х<1>~Х<1)жХ<2>~5жХа>+1 280 RETURN:END

Пример. В программе со строки 260 записывается подпрограмма вычисления правых частей (4.16) в виде

f(/)=(функция переменных Х(/)).

В данном случае подпрограмма служит для решения системы уравнений

x,+31gx,-xi=0, 2xi - jfijf2 - 5x1 +) = 0.

Задав N=2, M=10 (максимальное число итераций для случая, если решение расходится), 8=£=lIO- x,o=A:1(0)=3,4 и Х2о=Х2(0)=2,2, получим х, =3,487442788, Х2=2,261628631 и S=3. Значение переменной 5 дает число итераций, а значения Р{Г) - невязки системы. Время счета контрольного примера 16 с.

Решение (4.16) с помощью методов минимизации функций ряда переменных описано в § 4.6.

§ 4,5. Решение алгебраических уравнений с действительными и комплексными коэффициентами

Алгебраические (но ие трансцендентные) уравнения сводятся к канонической форме

где коэффициенты Ао, А\, An-в общем случае комплексные" числа. В таком виде алгебраические уравиеиия используются редко. На практике обычно встречаются полиномиальные уравнения с действительными коэффициентами

аоХ" + а,х- + ...+а„ ,х + а„=0, (4.18)

причем нередко (4.18) записывается в другой форме:

апХ" + ап-1х- + ...+а,х+ао=0. (4.19)

Отметим некоторые известные свойства алгебраических уравнений, на которых базируются методы их решения.

1. Корни алгебраических уравнений могут быть действительными и комплексно сопряженными, причем на комплексной плоскости все они лежат в пределах кольца, ограии-ченного окружностями с радиусами Л=1 + + амЕкс/ап и г=1 + (аСакс/ап1, где а„акс- наибольший по модулю коэффициент полинома, аяшс- следующий наибольший по модулю коэффициент (т. е. без учета а„акс)-

2. Число корней с учетом их кратности равно п, причем полиномы нечетной степени обязательно имеют хотя бы один действительный корень.

3. Число положительных (отрицательных) корней равно числу перемен (постоянств) знаков в последовательности а„, ап-\, ..... со коэффициентов или меньше его на четное число.

Квадратное уравнение

имеет решение в виде

Xl,2= -

a2X+aiX+ao=0

(4.20)

причем если подкоренное выражение отрицательно, корни получаются комплексно сопряженными, если оно положительно, оии действительные.

Программа 4.21.

Приме р. Для уравнения 2х-5х-10= = 0 получим х, = -1,311737691 и Х2=3,811737691, а для уравнения х + 2х + + 15 = 0 получим Х.2=-1±/-3,741657387. Приведенное кубическое уравнение x + a2x-f a,x-fao=0 (4.21)

может решаться по точным формулам [18]. Однако более компактное решение получается при использовании следующего алгоритма.

Вначале методом половинного деления с е=ЫО~° находится действительный корень Хз (он обязательно существует). Затем делением (4.21) на (х-хз) получается квадратное уравнение, из которого вычисляются два других корня xi и хг.

Программа 4.22.

Пример. Для уравнения х-6x+21х- -52=0 получаем х,.2=1 ±/-3,464101615 и хз=4, а для уравнения х-6х+11х-6=0 получим Xi=3, Х2 = 2 и хз=1.

Уравнение (4.19) с действительными коэффициентами решается методом Хичкока. При этом (вначале) итерационным методом из (4.19) выделяется квадратичный множитель (если rt2) и вычисляются два eto корня. Затем (4.19) делится на этот множитель, т. е. п уменьшается на 2, и процесс повторяется до вычисления всех корней.

Программа 4.23.

Пример. Для уравнения х+9х+ +31х*+59х+60=0 прлучим кории xi = = -1+/-2, Х2= -1-/-2, хз=-3+/-0, Х4=--4+/-0 (т. е. Хз и Х4 действительны), а для уравнения х+8х + 31х + 80х* + 94х + + 20 = 0 получим Х = -0,2679491924, Х2 = = -2, хз= -1+/-3, Х4= -1-у-З и Хб = = -3,732050808 (при вычислениях принималось 6=Ы0-*).

Решение уравиеиия общего вида (4.17) с комплексными коэффициентами Л=а/+/6,-обобщенным методом Ньютона реализуется следующим алгоритмом.

1. Вводим массив Л,- и задаем точность вычислений 6.

2. Задаем начальное приближение Zo= * -xo+jyo, в частности Хо=0 и уп=0 (Метод дает сходимость к какому-либо корню при любых начальных приближениях).

3. По методу Ньютона (ио е применением аппарата комплексной арифметики) находим корень (4.17) Zl.

4. Делим (4.1.7) иа {Z-Z,), т. е. понижаем степень (4.17) иа единицу и повторяем вычисления с п. 2 до нахождения второго корня и т. д. до тех пор, пока ие будут вычислены все кории.

Более подробное описание этого метода см. в. , [24]. При использовании програм-

Программа 4.22.

05 printрешение приведенного кубического уравнения

10 1НРиТвведите а2=с! шритвведите а1=в

15 inputвведите а0=а: inputвведите а rmV=k

20 letx=0sletk=k+1

25 letg=x+c!letl=6!«x+b

30 if L:«X+A>=0 then 50

40 letx=x+k

50 letk=k2!letx=x-k

60 if k>=1e-10 then 25

70 lete=-6/2!leth=e*e~l

80 if M>=0 then 130

90 printXb2=e+- JiKsqr(-M>

100 prihtx3=x!st0p

130 letf=SQRi:M) s leth=f+e: leti=e-f

140 printx1=h5 printx2=i:printx3=x!end

Программа 4.23.

05 PRINTВЫЧИСЛЕНИЕ ВСЕХ КОРНЕЙ ПОЛИНОМА 10 PRINTC ДЕЙСТВИТЕЛЬНЫМИ КОЭФФИЦИЕНТАЖ 15 ШРИТЗАДАЙТЕ ПОГРЕШНОСТЬ Е=Е

20 INPUTЗАДАЙТЕ СТЕПЕНЬ ПОЛИНОМА N=N!DIM А«:гжН+1>

25 FOR J=l ТО N+1SPRINTВВЕДИТЕ А!2.e!N+l-J

30 INPUT A<J>!NEXT J!LETK=1

35 LETT=1SLETC=A(2>/A<1)

40 IF N=1 THEN LETP=-C!LETG=0s6OTO 270

50 IF H=2 THEN LETH=C:*.C4-A<3)/A<1)!60T0 220

6ie LETri=l0sLETC=4! LETD=8s LETU=4

70 UETU=8! LETF=1: LETf)J=2i LETT=0

80 IF М<>10 THEH 100

85 LETP=C:LETM=0 s LETQ=B s LETC=U:LETD=U

90 LETU=PiLETU=Q!LETY=C!LETZ=D!LETF=-F

100 LETM=M+1!LETH=0! LETO=A( 1 >sLETP=A<2>-CiKQ: LETL=Q

120 FOR J=3 TO N!LETR=P!LETP=A<:J>~C*.R-D*e

130 t.ETQ=R!LETR=L!LETL=Q-CiKR-H*IisLETH=R:NEXT J

140 LETQ=A<N+l)-DiteGsLETS=L+CiiiR

150 IF T=0 THEN LETX=IiiKR:LETH=RiKX+SiKL!lF H=0 THEN 80

160 LETC=C+<PiKS-QiKR)/H!LETIi=D+<PiKX+QiKL)H

170 IF C-V+D-ZO0 THEN 190

180 IF Р=-Ы THEN PRINTНЕТ РЕШЕНИЯsSTOP

185 LETW=-F

190 LETH=C*C4-Iil IF S0R<<Q-P*.C2)2+PiiiPiiiABS<:H))>E/N THEN 200 LETT=0iLETA<2)=A<:2)-CiKA<l)!FOR J=3 TO N-l 210 LETft<J)=A<J)-C.iKA<J-l)-D*A<J-2)sNEXT J 220 LETP=-C/2!LET0=SQR(ABSCH))

230 IF H>=0 THEN LETM=P+Q!LETP=P-QsLETQ=0!6OTO 25© 240 LETM=P

250 PRINT!2.0!X<K)=!F1.9!M +JiK<Q>:LETK=K+1 270 PRINT!2.0!X<K>=!F1.9!P +JiK(-Q> !LETK=K+1 290 IF TO0 THEN 35© 300 LETN=N-2!LETA=0

310 IF SQR<<S-RiKC2)-2+RiKRiiiABSCH)><=E THEN LETA=1 320 LETB=0!IF N>=2 THEN LETB=1 330 IF A+B=2 THEN LETT=l!60T0 110 340 60T0 35

350 PRINTВЫЧИСЛЕНИЯ ЗАКОНЧЕНЫsENE

05 PRIWTРЕШЕНИЕ КВАДРАТНОГО УРАВНЕНИЯ

10 ШРИТВВЕДИТЕ А2=С! ШРиТВВЕДИТЕ А1=В

15 ШРИТВВЕДИТЕ А0=А

20 LETri=B-2/C5LETF=DiKri-A-C,

30 IF F>0 THEN ,60

40 PRINTXb2=-ri4- JxSQR(.-f>

50 60T0 10

60 PRINTXl=~ri~SQRi:F> . 70 PRINTX2=-D+SQR(F> 80 GOTO 10:ENB