Пример. Для системы дифференциальных уравнений и данных к программе 4.53 будем получать:

0,05

0,15

0,25

2,499985857-10" 9,999968178-10-- 2,249994626-10" = 3,999991927-10- = 6,249988623-10- =

2,049999950 2,099999890 2,149999819 2,199999736 2,249999637

Время счета для одной точки около 10 с, т. е. вдвое больше, чем у программы 4.53. Однако вдвое больше и выбранный программой шаг h В целом скорость счета у программ 4.53 и 4.54 примерно одинакова и зависит от характера решаемой системы дифференциальных уравнений.

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

Многошаговые методы решения дифференциальных уравнений базируются на использовании данных реше51ия на нескольких предшествующих шагах. Это позволяет по-

выйть скорость вычислений. Однако для начала вычислений приходится выполнять одношаговыми методами несколько первых шагов. Аналогично это делается при каждой смене шага интегрирования.

Методы прогноза и коррекцшг - разновидности многошаговых методов, при которых решение вначале проводится по формуле прогноза, а затем уточняется [ю формуле коррекции.

Ввиду сложности программной реализации многошаговых методов, а также неявных методов численного интегрирования (последние лишены численной неустойчивости решения), они очень редко используются при решении задач на персональных ЭВМ. .Алгоритмы их реализации описаны в [10, 18, 24, 36].

§ 4.11. Гармонический синтез

гармоничесшй синтез - получение колебаний сложной формы путем суммирования их гармонических составляющих (гармоник). Гармонический синтез может иметь как самостоятельное значение, так и использоваться как средство тригонометрической интерполяции функций (сонместно со спектральным анализом - см. § 5.1).

65 PRINTРЕШЕНИЕ СИСТЕМЫ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ

10 PRINTМЕТОДОМ РУНГЕ-КУТТА-ФЕЛЬБЕРГА С АВТОМАТИЧЕСКИМ

15 PRINT ВЫБОРОМ ШАГАsPRINT!2.0!

40 INPUTВВЕДИТЕ ЧИСЛО УРАВНЕНИЙ N=N

45 DIM V<N>/L<:N),W<N)rA<:N)fB<N>fCaO;riaOfEaOrF<:N)

50 INPUTЗАДАЙТЕ ПОГРЕШНОСТЬ ВЫЧИСЛЕНИЙ Е=Е1

60 INPUTЗАДАЙТЕ НАЧАЛЬНЫЙ ШАГ Н=Н

70 INPUTЗАДАЙТЕ НАЧАЛЬНОЕ Х0=Х

80 FOR J=l ТО N SPRINTВВЕДИТЕ НАЧАЛЬНОЕ V0<:J)=

96 INPUT w<j)!letv<:j)=u<:j)snext j :

110 LETE3=0s6OSUB 4Ci0!LETri=0sFOR J=l TO N

120 LETA<:J)=F<J»ttH!LETV<:j>=W<J>+2)ttA<:J)/9!NEXT J

130 LETX=X+2)t!H-9!60SUB 460! FOR J=l TO N

140 LETB i J > =НжЕ С J ) s LET V С J ) =Ы С J ) +A С J > .- 1 £+B С J > ,4: NEXT J

150 LETX=X+H/9!80SUB 406: foi? J=l TO H

160 LETC«:j)=H)ttF<:j)!LETV<:j)=w<:j)+<34.5«A<:.j)-i2i.5*B<:j>+i35-*c<:.j)>-64.

170 NEXT .JsLETX=X+H/2.4:60SUB 406: FOR J=l TO N

180 LETD С J >=НжЕ С J ) s LETG=W С J ) -17Ж A С J >12+6.75*B С J >

185 LETV<:J>=G-5.4»C<:j>+16)ttri<:J>-15:NEXT J

190 LETX=X+H/4!60SUB 400:for J=l TO N

195 LETE<:a)=H»F<J>sLETG=W<:J)+65«A<:J)---432-.3125!«Bf,J>

200 LETV<:J)=e+.8125«C<:J>+4»ri<:J)/27+5«E<:j>-144!NEXT J

266 LETX=X-H/6s60SUB 400!FOR J=l TO N

210 LETl. С J ) =H»F С J ) s LETV С J ) =W С J ) +A ( J > /9+. 45«C < J )+1 €m < J ) -45+E С J ) 12

220 LETG=A<:J)/150-.03»C<J)+16«D<:J>-75+E<:.J>-£0

230 LETE2=ABS<Q-.24»L<:J>):IF E2<=E1 THEN 250

240 LETE3=l!60T0 260

250 IF E2<El/20 THEN LETri=D+l

260 NEXT JsLETX=X+H/6s60SUB 400:IF E3=0 THEN 296 270 LETX=X-H!FOR J=l TO NsLETVCJ>=W<J>!NEXT J

280 LETH=H/2!60T0 118 . -

290 IF ri=N THEN LETH=H+H

3C10 РР1НТ!Е1.9!ДЛЯ X=X!for J=l TO N

310 PRINT!2.0!V<: J>= !F1.9!V<:j>

320 letw<:j>=v<:j)!NExt J!6oto ii6

400 LETFa)=V<l>+V<:2)-X)ttX+X-2

410 LETF<:2)=-2)ttVa>+4)ttV<:2>+2)ttX)ttX-4«X-7

420 fTURNsENn

Гармонический синтез с помощью ряда Фурье (1 вариант) заключается в вычислении ограниченного т гармониками ряда Фурье

У (О ="2-+ X 2"?! + * -" 27ikf,t),

* = i

где t - время, f\ - частота первой гармоники, А/, к Bk - коэффициенты ряда Фурье. Программа 4.55.

§ 4.12. Вычисление собственных - значений и векторов матриц

Собственные значения К (t = 1, i, ..., п) квадратной матрицы А = [а,,] п-го порядка есть числа (действительные или комплексные), удовлетворяющие условиям Лх = Х,х. В механике собственные значения характеризуют, например, внутренние напряжения в объекте, находящимся под воздействием

10 PRINTГАРМОНИЧЕСКИЙ СИНТЕЗ 1

20 INF4JTВВЕДИТЕ ЧИСЛО ГАРМОНИК M=M:DIM Й<ЮгВ<Ю

30 РКШТВВЕДИТЕ А 6: INPUT А

40 FOR N=1 ТО M:PRINT!2.0!ВВЕДИТЕ AN

50 INPUT АаОSPRINTВВЕДИТЕ BNsINPUT BCNsHEXT Н

60 INPUTBBEAHTE ЧАСТОТУ F=F

70 INPUTBBEAHTE ВРЕМЯ T=TsLETC=2»*PI»F)ttT!LETS=0

ее FOR N=1 TO M: l.£TD=C»N

90 LETS=S+<A<N)»C0S<D>+B<N)»SIN<:iO)sNEXT N

100 LETV=A/2+SsPRINTV<T)=!Fl,9!V!e0T0 rosEND

Пример. Для Ло = 4, Л1 = 10, Bi=5, /42 = 8, В2 = 2, A:i=4, Bs=l. 4=1, В4=0,2, 5 = 0,1, В5 = 0,05 и fi = l будем получать: у (0)=25,1; у (0,1) = 16.32687379; у (0,2) = = 0,9348395361.

Гармонический синтез с помощью ряда Фурье (2 вариант) сводится к вычислению ряда

У (0=-т+ У А, COS (2nkf,t + Qk),

внешних усилий, действующих со всех сторон. В электро-радиотехнических устройствах собственные значения матриц определяют характеристические постоянные времени и режимы работы этих устройств. Все это говорит о практической важности задач на нахождение собственных значений (4, 24, 41].

Определитель

D (Х) =

заданного гармониками с амплитудрй А/, и фазовым углом (рад). Программа 4.56.

01 -X а\2 ... а,„

0.21 022 -- - . . 02„

. .. о„„ - X

называется характеристическим многочленом

10 PRINTГАРМОНИЧЕСКИЙ СИНТЕЗ 2

:20 INPUTBBEAHTE ЧИСЛО ГАР110НИК M=MsriIM AOUfQan)

30 PRINTBBEAHTE А 0sINPUT A

40 FOR N=1 TO M!PRINT!2.0!ВВЕДИТЕ AN

.50 INPUT A(:N>:PRINT!2.0!ВВЕДИТЕ ФАЗУ В РАДИАНАХ QN

55 INPUT eCHisNEXT N

6Й INPUTBBEAHTE ЧАСТОТУ F=F

70 INPUTВВЕДИТЕ ЕРЕМЯ T=T!LETC=2»#PI»F»T!LETS=0 80 FOR ft=l TO MsLETS=S+Aao«COSfC»N+Q<N>>sHEXT N 100 LETV=ft£+S! PR INT V С T >= ! F1.9! V! БОТО 70 s END

Пример. Для Лс = 4, Ai = \0, Qi=0, 2 = 8, Q2 = 0,5, Аз = 5, (Зз = 0,8, У!4 = 2, Q4 = 2, 5 = 0,5, Q5 = 3 и /, = 1 будем получать: у (0) =21,17690412; у (0,1) =4,223677837; у (0,2) = -2,576641635.

Аналогичным образом осуществляется гармонический синтез при других тригонометрических рядах. Например, при ряде Котель-никова с ограниченным числом членов [18]

(/) = £ y(kAt)

k = 0

sin л (t/At - k) п (t/M - k)

где At - шаг, с которым задаются отсчеты у (k At). В некоторых случаях гармонический синтез упрощается. Так, если предполагаемая функция у {{) симметрична относительно оси /, то амплитуды всех четных гармоник равны нулю.

(полиномом) матрий,ы А. Уравнение Ах = %х эквивалентно уравнению (А - 1.Е) х=0. Уравнение (1е1(Л-Л£)=0 называется характеристическим многочленом матрицы А:

det (/4-Х£)=Р (Х) =

= -о,-02Я, -о.зХ- ... -о„Г-+Г = 0.

Его корни определяют все собственные значения Х;.

Каждому h соответствует принадлежащий ему собственный вектор Yi с компонентами У12, Ут- Они удовлетворяют решению системы уравнений AYi = hYi.

Коэффициенты о, характеристического полинома Р (Л) находятся методом Данилевского. Для этого рассматривается специальным образом преобразованная последова-

тельность матриц АК А-К

при-

чем преобразования (см. ниже) направлены на то. чтобы получить матрицу А , у которой коэффициенты а,- характеристического полинома оказываются ее последним столбцом. Преобразования Л* в Л"""" выполняются по формулам

где ВП>.1АП<.!/АП>л и B(f>=(«-

-A?M%,.i {1фк+\), причем /=1, 2-----N;

/=1, 2.....N. Для экономии памяти ПЭВМ

вспомогательная матрица В заносится на место матрицы А. Программа 4.57.

мые программрй. Степень первого полинома равна первому значению k, степень второго - разности второго и первого значений k и т. Л (степень последнего полинома равна разности N = 11 и последнего значения k). Так, если программа выдает k, = 2, fe-. = 3 {N = b), а массив а, выдается в виде пяти чисел

Ш щ as

I \ I

dS akV cQ,

то степень первого поЛинома k\ = довательно,

Р, ()==-<-4У?. + ;.=.

= 2 и, сле-

05 PRINTВЫЧИСЛЕНИЕ КОЭФФИЦИЕНТОВ ХАРАКТЕРИСТИЧЕСКОГО ПОЛИНОМА 10 PRINT МАТРИЦЫ МЕТОДОМ ДАНИЛЕВСКОГО

20 1НРиТЗАДАйТЕ ПОРЯДОК МАТРИЦЫ N=N:D1M ft(N,H> 25 FOR 1=1 ТО NSFOR J»l TO N

30 PRIHT!2.0!ВВЕДИТЕ AI/JsINPUT A<bJ>sHEXT JsNEXT IsLETQ=l

40 FOR K=Q TO N-l«LETL=0sFOR I=K+1 TO N

50 IF A8S<L>-ABS<A<bK))<0 THEN LETL=A< bK>!LETP=I

60 NEXT IS IF LO0 THEN 90

70 FOR J=e TO KILETA<0/J>=A<J,K)SNEXT J

80 LETQ=K+1IPRINT!2.0IK=KS6OTO 40

90 IF K+1=P THEN 140

100 FOR J=K TO NsLETR=A<K+bJ>

110 LETA<K+bJ)»A<P/J)sLETA<P,J)=RsNEXT J

120 FOR J=Q to N:LETR=ft<J,K+l>

130 LETA<J,K+1>=A<J/P)SLETA<J,P>=RSNEXT J

140 FOR J=Q TO NsLETA<0/J)=A<J/K)sNEXT J

150 LETC=A<K*1/K)SF0R J=K TO M

160 LETA<K+l/J>»A<K+l/J)/CsFOR I=Q TO N

170 IF I=K+1 THEN 190

180 LETA<l,J)»A<bJ>-ft<0,l)*ft<K+l/J>

190 NEXT UNEXT J

200 FOR 1=Q TO NsLETS=0«FOR J=K+1 TO N

2ie LETS«S+A<I,J)!«A<0/J)sHEXT J

Z№ IF I-K-K THEN 240

230 LETI>=0:6OTO 260

240 IF I=Q THEN LETD=0s6OTO 260

250 LETD«A<0,I-1>

260 LETA<bK+l)=S+D«NEXT IsNEXT К

270 FOR J=e TO NsLETA<0/J>=A<J/N>sNEXT J

280 FOR. J»l TO NlPRINT!2.0!AJ!F1.9! = A(0/J)

290 NEXT JSEND

Пример. Если A - невыро.жденная матрица, программа 3.69 дает значения Л/ = 0 характеристического полинома. Например, для матрицы

-5,509882 1,870086 0,422908 0,287865 -11,811654 5,7119 0,(149099 4,308033 -12,970687 L 0,006235 0,269851 1,397369

Степень второго полинома ke - k\ = \ и

Р-2 il) = -a4i+%. Степень третьего полинома N - k-2 = 2 и

Для матрицы

получаем /11 =0, = -12296,55057, А2 = т = = -5349,455515, /13=а, =-797,2787648 и /14=04=-47,88843.

Если А - вырожденная матрица, то массив Ai = ai содержит коэффициенты с обратными знаками всех составных полиномов. Идентификатором полиномов в этом случае являются значения переменной /г, выдавае-

программа 4.57 выдает k = k\=2, Л 1=0, Л2 = 2, Л3 = 2 и Л4 = 5. Составляя полиномы по описанным для этого случая правилам, получим

/5, (Я.) = О - 2>.= - 2Х + X-

Pi (>,) = -2 -5>L-fX=.