(Q,Z) (1,Z) (Z,2)

(2,1)

a,o)

В as

Phc. 4.11. Расположение узлов при вычислении двойных определенных интегралов методом Симпсона (а)

и Гаусса (б)

несколько квадратов, вычислив /, для каждого j-ro квадрата. .Программа 4.46 предусматривает последовательное суммирование значений /,-. Например, разбив область интегрирования (4.31) на 4 квадрата (а==4; 6,=4,2; с,=2,3; d,=2,6; 02 = 4,2; б2 = 4,4; С2=2,3; rf2 = 2,6; аз=4,2; 6з = 4,4; Сз = 2; йз=»2,3; 04 = 4; Ь, = 4Х « = 2 и 4 = 2,3 - см. рис. 4.11, б), получим

/=0,02500605059.

Время счета для каждого квадрата менее 2 с.

Кубатурная формула Гаусса с п = 2 имеет вид

/=5 S f {X, y) = PR и (X,, {/,) +f (Х2, {/2) + о с

+ Мх.з,4/з)+/(-«4,1/4)1, (4.32)-

где Я=(6-а)/2, R={d-c)/2, а узлы интерполяции расположены на расстоянии ±КР и +KR от осей, где K=l/.

Программа 4.47.

Пример. Для интеграла (4.31) с областью интегрирования в виде одного квадрата получим / = 0,02500532089, а с областью интегрирования в виде четырех прямоугольников - / = 0,02500594169. Время вычисления для одного квадрата области интегрирования - около 1,5 с.

Сложная кубатурная формула Г аусса реализуется разбивкой отрезков [о, Ь] на М "и [с, d] на N частей с применением для каждого элементарного прямоугольника формулы (4.32).

Программа 4.48.

Пример. Для интеграла (4.31), задав M=/V = 4 (всего 16 прямоугольников), получим / = 0,02500598247 при /с = 24 с.

Аналогичным образом могут быть составлены программы для вычисления интегралов с большей кратностью. Формулы численного интегрирования для них приведены в [18, 36]. Интегралы большой кратности вычисляются методом Монте-Карло (см. [4]).

Программа 4.47.

18 PRINTВЫЧИСЛЕНИЕ ДВОЙНОГО ИНТЕГРЙЛЙ ПО ПРОСТОЙ

15 PRINT ФОРМУЛЕ ГйУССЙ:ЬЕТ8=е

20 INPUTВВЕДИТЕ ПРЕДЕЛЫ Й/В/С/Г Й,В,С,Ъ

36 ЬЕТ0=<й+Е)/2! LETP=<B-ft)/2! LETe=<C-H)>/2

40 LETR=<D-C)/2!LETK=SGiRa/3)

50 LETX=0-bP*K!LETV=Q-i-R*K!60SUB 168

60 LETX=0-P*K!60SUB 166

70 LETV=Gi-R*K:&OSUB 166

88 LETX=0-frP*K!60SUB 166

90 PRINTЗНАЧЕНИЕ ИНТЕГРЙЛЙ 1=8жРжр!&0Т0 26 180 LETF=l>X/V!LETS=S-bF:RETURN!ENIi

Программа 4.48.

18 PRINTВЫЧИСЛЕНИЕ ДВОЙНОГО ИНТЕГРАЛА ПО СЛОМНОй

15 PRINT ФОРМУЛЕ ГАУССА!LETS=e

20 INPUTВВЕДИТЕ ПРЕДЕЛЫ A/B/C/D Ы,В,С,Ъ

25 WUTВВЕДИТЕ ЧИСЛА Mj-N M,N

30 LETH=<B-W)/M: LETL=<Ii-C)/N! LETK=SGiR< 13)

40 FOR 1=1 TO N! LETIi=C-bL

50 LETe=<C-Hi)/2: LETR=L/2!LETA=W

60 FOR J=l TO M:LETB=A-bH:LET0=<A-KB)/2!LETP=H/2

70 LETX=0-bP*K:LETV=Q-i-R*K:60SUB 140

80 LETX=0-P*K:&OSUB 140

90 LETV-Q-RsifK:60SUB 140

100 LETX=0-bP»K!60SUB 140

110 LETA=BSNEXT J

120 LETC=Ii:NEXT I

130 PRINT3HA4EHHE ИНТЕГРАЛА 1=S*H«L/4:ST0P 140 LETF=1XV! LETS=S-bF: RETURN! END

§ 4.10. Решение систем дифференциальных уравнений

Задача Коши заключается в решении систем обыкновенных дифференциальных уравнений первого порядка, представляемых в виде dy,

-F\ (х, у....., у,.....у).

= Ft (х, у,.....у,.....{/д,), (4.33).

dyN dx

= Ff, (х, yj.....j/д,).

где / = l-e-yv - номер каждой зависимой переменной у,, х - независимая переменная. Если задача Коши решается для анализа поведения системы или объекта во времени, то X является временем (х = /). Решение системы (4.33) при заданных начальных

условиях Х = Хо, yi (Хо) =У\0, г/2 (Хо) =i/20, ..

.... Ук(хч)=Ут сводится к нахождению зависимостей (интегральных кривых) yi (х), t/2(x), У; (х), </yvW. проходящих через точки, заданные начальными условиями (хо, у\о), (хо, У20)..... (хо, У10) ,

(хо, {/дго)- Задача Коши сводится к интегрированию дифференциальных уравнений. Порядок метода численного интегрирования при этом определяет и порядок . метода решения (4.33).

Обобщенная форма записи каждого из уравнений системы (4.33) может быть представлена в общем виде

=f,(x.K.).

где Y; в правой части уравнения - векторы

переменных у>, уг.....у,.....у, а f,•-

правая часть каждого из уравнений (4.33). В частности, одно дифференциальное уравнение (j/=Kj = Ki и Fj-F=F,) записывается в виде

.я.

Дифференциальные уравнения высшего порядка

= f (X, у, у, у".....у"- "), (4.35)

где (п) - порядок уравнения, могут, быть сведены к системам вида (4.33) или . (4.34) с помощью следующих преобразований:

dyn-2 dx

=F(x.y.y......</„-,).

(4.36)

Следовательно, решение (4.35) сводится к решению системы дифференциальных уравнений первого порядка (4.36).

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

y,(,-Hn = yj,-f/if/ (х„ К,-.),

где h - шаг интегрирования (приращение переменной х). Этот метод обладает большой погрешностью и имеет систематическое накопление ошибок. Погрешность метода R ~ (fi), т. е. пропорциональна И..

Метод Эйлера - Коши с итерациями заключается в вычислении на каждом шаге начального значения

jioVo=Ki, + Afi (л, К;,). Затем с помощью итерационной формулы У<г+о=У;/+- [Fi (Xi, Yji)+F; (х,+,, Yt+\>A

решение уточняется. Итерации проводят до тех пор, пока не совпадает заданное число цифр результата на двух последних шагах итераций. Погрешность метода R-(h). Обычно число итераций не должно превышать 3-4, иначе нужно уменьшить шаг ft.

Модифицированный метод Эйлера второго порядка реализуется следующими рекуррентными формулами:

где Yfi+,/2)=Yji + hFj (xi, Yi,)/2. Метод дает погрешность R -~ (А) и имеет меньшее время вычислений, поскольку вместо нескольких итераций производится вычисление только одного значения Kf(, + i/2)-

Метод трапеций - одна из модификаций метода Эйлера второго порядка. Он реализуется применением на каждом шаге формулы

где Ki\hFi(Xi, К,,), Kii=hF, (Xi+h. Yji+ + Ki\), и дает погрешность R~(h). Этот метод относится к общим методам Рунге - Кутта (см. далее). Программа 4.49.

Пример. Вычисление функций F, (х, Yj) оформляется подпрограммой, записываемой со строки 200 и завершаемой операторами RETURN и END. При этом форма записи каждой функции следующая:

F (J) = {функция X,Y (1), Y (2).....Y (N)>,

где J - номер, соответствующий номеру / в обобщенной форме записи Fj (х, yj). На-. пример, для решения системы дифференциальных уравнений dy.

dx dy2 dx

-=f, (х,уиу2)=У1,

= F2 (X, уи У2) =(--.2) -]р-У>

подпрограмму следует записать в виде F(1)=Y(2),

F(2) = (Y(1)/X-Y(2))/X-Y(1)

(см. также выражения в строке 200 программы 4.49). Для N=2; /г = 0,1; хо = 0,2; 4/10=0,099500833 и I2o = 0,49235 получим:

X у,

0,3 0,148367183

0,4 0,1961203205

0,5 0,2424068609

0,440372468

0,4830282264 0,4701529365 0,4537779586

0,3249711437

и т. д.

Время счета на каждом шаге около 1,5 с. Не рекомендуется применять этот метод, если возмо.жны изменения у, (х) с сильно различающейся крутизной (например, одни

/iC3;=ftf/(x,-f,K,-,-fi К2,); K4i=KFi(xi+h, К„ + /Сз,); К,(,+п= >,;+ -i- {Ku+2K2i+2Ksi+K4i).

При переходе от одной формулы к другой задаются или вычисляются соответствующие значения х и К; и находятся по подпрограмме значения функций F, (х. У,).

Решение одного дифференциального уравнения методом Рунге - Кутта производится по приведенным формулам, если в них опустить индекс /, а из алгоритма исключить циклы.. Последнее резко упрощает программу и позволяет получить минимально возможное время счета.

Программа 4.50.

10 PRINTРЕШЕНИЕ ОДНОГО ДИФФЕРЕНЦИАЛЬНОГО УРАВНЕНИЯ 20 PRINT МЕТОДОМ РУНГЕ-КУТТА

30 INPUTЗАДАЙТЕ ШАГ Н=Н!INPUTЗАДАЙТЕ НАЧАЛЬНОЕ Х0=Х

40 INPUTВВЕДИТЕ НАЧАЛЬНОЕ V0=Z:!LETV=Z:

50 60S;UB 100sLETA=F«HsLETK=X+H/2!LETV=Z+A-£

60 60SUB 100sLETB=F«H5LETV=Z+B/2

70 &0SUB 100!LETC=F«H!LET!«;=X+H/2:LETV=Z+C.

80 60£Ш 100sLETV=Z+(A+2«<B+C.)+Hs*F>/65LETZ=V

90 РР1НТДЛЯ X=X5PRINTV=Vs60T0 58

100 LETF=-V! RETURN! ENII

участки кривых у, (x) плавные, a другие резкие), поскольку, при этом возможно возникновение численной неустойчивости решения. Ее устранение может потребовать выбора очень малого шага h, что значительно увеличивает общее время вычислений.

Метой Рунге - KyiTa четвертого порядка является наиболее распространенным методом решения систем (4.34) при шаге h = = const. Его достоинством является высокая точность - погрешность R -~ (h) - и меньшая склонность к возникновению неустойчивости решения. Алгоритм реализации метода Рунге - Кутта заключается в циклических вычислениях Уц1+\\ на каждом г+1 шаге по следующим формулам:

K,i = hFi(xi, K„);K2; = Af,(x.+g , У„+1 K,j);

Пример. Для дифференциального уравнения у= -у/т= -у при т=1, /г = 0,1, хо = 0 и 4/0=1 будем получать у (0,1) = = 0,9048375, у (0,2) =0,8187309014, у (0,3) = = 0,740818422, у (0,4) =0,6703202889"

у (0,5) =0,6065309344. Время выдачи одного значения у (х) около 1 с. Подпрограмма вычисления

F=.dy/dx

записывается со строки 100.

Реализация решения системы из N дифференциальных уравнений иллюстрируется программой 4.51. Максимальное число зависит от сложности уравнений и объема ОЗУ ЭВМ. Для ПЭВМ типовые максимальные значения N лежат от 10-20 (для ПЭВМ класса Pocket Computers) до 50-100 (для профессиональных ЭВМ).

le PRIHT-РЕШЕНИЕ СИСТЕММ ДИФФЕРЕНЦИйЯЬНМХ

20 PRINT УРАВНЕНИЙ МЕТОДОЙ ТРАПЕЦИЙ

30 INPUTЗАДАЙТЕ ЧНСПО УРАВНЕНИЙ N=N

40 liiri V<N>,W<N>rKaO.-F<N>

100 INPUT ВВЕДИТЕ НАЧАЛЬНОЕ Х0=Х

110 FOR J=l ТО NSPRINT!2.0!ВВЕДИТЕ НАЧАЛЬНОЕ VOCJJ

120 INPUT W.<J):LETV<:J>=UKJ>sNEKT J

130 INPUTBBEAHTE ШАГ H=H

140 60SUB 200:FOR J=i TO N

150 LETK;<.J)=Hi«F<J>-.LETV<J>=W<.J>+K<J>!hEXT J

160 LETX=!«;+HsPRINT!F1.9!flnS X=Xs60SUB 280

170 FOR J=l TO NsLETVf,J)=W<J>+<:K<:.J>+H«F<:J>>2

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

190 LETW<J>=V<J>!NE!«;T JS60T0 140 ,

200 LETF<l>=V<2)sLETF<:2>=<V<l)K-V<:2>)/X-Va)

210 RETiJRNsENn