а затем для расчета интеграла (2.20) суммируют Ij по всем интерва

Перейдем к рассмотрению наиболее употребительных кватп ных формул. Ц.

Формула трапеций. При ее получении функция / (х) ннтепп руется на каждом элементарном отрезке [xj, Xj+,] линейной функ,?" (рнс. 2.10). Тогда -f

ifiXj)-f{XHl)] {ХШ~Х}1

Xi =a, Xn+i =6. (2.25

и после суммирования по всем интервалам ixj, Xj+] придем к форь

Формула Симпсоиа. Она основана на использовании отрезка разбиения с тремя узлами и интерполяции подынтегральной функи™ полиномом второй степени по трем равноотстоящим узлам а,

fix)

Рис. 2.JO

Рис. 2.11

Xj+2- На отрезке Uj, Xj.+ 2h] функция f (х) заменяется параболой, проходящей через точки fj, fj+i, fj+ (рис. 2.11). Интеграл от интерполяционного полинома на отрезке [xj, Xj + 2h] равен

Л= J fix)

Считая, что полное число интервалов N = (Ь - a)lh четное, и сУ мируя по всем парам интервалов, получим формулу Симпсона

(2.2

N/2 N/2

Л/3.

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

Теперь рассмотрим случай, когда узлы Xf не фиксированы, и та-м образом в квадратурной формуле (2.22) можно выбирать не толь- весовые коэффициенты с, но и расположение узлов Xi, в которых вычисляется подынтегральная функция. Использование этих дополнительных «степеней свободы» позволяет повысить точность квадратурных формул.

формулы Гаусса. Ьсли на отрезке интегрирования в качестве искомых параметров квадратурной формулы рассматривать (Л + 1) коэффициент Ci и (Л + 1} узел л, то получим (2Л + 2} неизвестных. Эти параметры можно выбрать так, чтобы квадратурная формула была точна для любого многочлена степени не выше (2N + !) Решение такой задачи известно: расположение узлов вычисляется с помощью корней полиномов Лежандра. Узлы Xi и весовые коэффициенты для различных приведены в [21. Построенные таким образом квадратурные формулы называют формулами Гаусса или формулами наивысшей алгебраической точности. Для гладких функций эти формулы дают очень высокую точность.

Оценка погрешности численного интегрирования. Различают два вида оценок: априорные и апостериорные. Априорную оценку получают заранее, до проведения расчетов, на основе теоретического анализа квадратурной формулы. Апостериорную оценку определяют после вычислений на основе сопоставления результатов расчетов, проведенных при разных числах отрезков разбиения.

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

R!, = U{x)6x-jI Cifixi) (=1

для которых удается найти верхние границы, справедливые при/i-»-(Л-к оо).

Определим такие границы для погрешностей формул прямоугольников и трапеций, используя разложение функций f (х) на отрезке XfJ в ряд Тейлора около точки х,- = (Xi + л;+])/2 и ограничи-аясь членами второго порядка:

f(x} = fCxi) + f4xinx-Xi)+rCx{x-~Xt)l2-i-.... Тогда для формулы прямоугольников получим

гр f{x)dx~[{xt)hP (~Xi]hV24,

а для формулы трапеций

т. е. погрешность вычисления интеграла иа одном элементарном нн тервале [х;, Xj+il пропорциональна Я.

Для определения погрешности вычисления интеграла на отрезк [й, Ь\ следует провести суммирование по всем подынтервалам:

-h2(6--a)fi,U/24. где fax - max [f (01. Аналогично для формулы трапеций получим

Показатель степени шага h в оценке погрешности называют мо. рядком точности квадратурной формулы. Квадратурные форму,!1ы прямоугольников и трапеций имеют второй порядок точности для функций, имеющих вторую производную, т. е. R -0 ф?). Заметим, что порядок точности обеих формул одинаков, хотя в одной использована интерполяция линейными функциями, а в другой - кусочно-постоянными.

Можно получить аналогичным образом выражение для главного члена погреишостн формулы Симпсона

> max тг,у

niax

т. е. формула Симпсона имеет четвертый порядок точности и является точной для полиномов степени не выше третьей.

Приведенные оценки неудобны для практического анализа погрешности, во-первых, потому, что определение значений f]

fmax может предстзвлять собой самостоятельную достаточно громоздкую задачу, а во-вторь[Х, эти оценки могут дать сильно завышенную величину погрешности.

Рассмотрим основной, применяемый на практике, способ апостериорной оценки погрешности, называемый правилом Рунге. Пусть из теоретического анализа известно, что численный метод имеет порядок точности р, т. е. погрешность R пропорциональна h:

и предположим, что постоянная А не изменяется при изменении величины шага h.

«2

max

Рсли получены значения Р и 1 в результате расчетов с шагами fi и 2 можно записать равенства:

I~nAhP, lPAhl.

Исключая неизвестное точное решение /, получим

A = {P-P)l{h\-h),

н тогда погрешность расчета, выполненного с шагом h.2, можно оценить по формуле

/ Я (Г- I)}\[h}hf ~ 11. (2.25)

Если бы оценка (2.25) была точной, то можно было бы найти и точное значение интеграла /. Однако изложенный способ обычно позволяет лишь правильно найти порядок величины погре[иности. Хотя иногда идея комбинации решений, полученных на разных сетках, и используется для уточнения решения.

Наиболее часто вьшолняют расчеты с удвоением шага, тогда hi 2К и

/:(Я„/1)/(2Р -1). (2.26)

Идею правила Рунге можно применять также для получения оценок погрешностей решений дифференциальных уравнений. В частности, на ее основе выводится приведенная в главе 1 формула (1.60) для полной погрешности численного решения обыкновенного дифференциального уравнения, в которой используются два численных решения, полученные на сетках разной густоты. При решении многих сложных задач такой путь оценки погрешности численного решения - единственно возможный.

Программное обеспечение для вычисления интегралов. Для численного интегрирования имеется достаточно обширное программное обеспечение. Разумеется, для того, чтобы реализовать вычисления по формуле прямоугольников (2.21) или по формуле Симпсона (2.24) с заданным шагом h, нет необходимости в поиске соответствующей стандартной подпрограммы, так как их нетрудно запрограммировать и самому. Ниже в качестве учебного примера (рис. 2.12) приводится подпрограмма вычисления интеграла по методу Симпсона, <оторая будет использована далее при реализации аналитического Решения (2.13). Необходимые для ее понимания сведения даны в ком-нтариях к тексту. Отметп.м лишь, что среди формальных параметров подпрограммы присутствует имя подпрограммы - функции задающей подыинтегральное выражение, поэтому в вь[-?,ьшающей программе это имя должно быть описано в операторе

.ternal.

некоторых случаях вопрос выбора квадратурной формулы и

1 с ПРОГРАММА РАСЧЕТА ОДНОМЕРНОГО НЕСГАЦИОНАРНОга ТЕИПЕРАТтюГО

2 С ПОЛЯ ПЛАСТИНУ ПО ТОЧНОМУ РЕЖНИЮ

3 С ИСХОДНЫЕ ДАННЫЕ :

4 С DL-ТОЛШИНА, АШ-ТЕПЛОПРОВОДНОСТЬ, АТ-ТЕНПЕРАТУРОПРОВОДНОСТЪ

5 С ALF-КОЭИИШЕНТ ТЕПЛООТДАЧИ, Т0-НАЧАЛЫШЙ ПЕРЕГРЕВ

6 С ПОДПРОГРАММА-ФУНКЦИЯ OV(X) ЗАДАСТ РАСПРВДБЛЕНИЕ МОЩНОСТИ

7 С Ш-ЧИСЛО ЧЛЕНОВ РЯДА

в С М-ЧИСЛО ИНТЕРВАЛОВ ПРИ ВЫСШйй ИНТЕГРАЛА

9 С П-ЧИСЛО ПРОСТРАНСТВЕННЫХ ТОЧЕК, JJ-ЧИСЛО МОМЕНТОВ БРЕЗйЩ

10 С Х(П)-МАССИВ КООРДИНАТ ПРОСТРАНСТВЕННЫХ ТОЧЕК

11 С TAU(JJ)-HACCHB МОМЕНТОВ ВРЕМЕНИ

12 DIMENSION A{50).BI{50),B2(50),C{50).T(10).X(10).TAU(2e) J3 С ОПИСАНИЕ ВНЕШНЕЙ ФУНКЦИИ W1NT-

14 С -ПОДЫНТЕГРАЛЬНОЙ ФУНКЦИИ ИЗОБРАЖЕНИЯ РАСПРЕДЕЛЕНИЯ НЙЦНОСШ

15 EXTERNAL iINT

18 С ОБЩИЯ БЛОК ДЛЯ ПЕРЕДАЧИ ПАРАМЕТРОВ В 17 С ПОДПРОГРАНМУ-ФУНКЦИС WINT

IB COMMON /INT/ AN.DL

19 С ВВОД ИСХОДНЫХ ДАННЫХ

20 READ I.DL,ALAM.AT,ALF,Ta

21 I FORMAT(5FI0.3)

22 PPJNT 2,DL,ALAM,AT.ALF,T«i

23 2 FORMAT{ DL=,G10.3, AUM=.С1в.З, АТ-.С1в.З/

24 ALF=.C10.3, T0=,F6.2)

25 REAI 3,NN,M,n,JJ

26 3 F0RMAT(4I5)

27 PRINT 4,fffl,M

28 4 FORMATC NN=,13, M=.I3) 2S READ I,(Xd),1=1,II)

30 PRINT 5,(X(l).M.in

31 5 FORMATC КООРДИНАТЫ ПРОСТРАНСТВЕННЫХ ТОЧЕК г7(5С11.3))(

32 READ I,(TAU(J).J=1,JJ)

33 PRINT 8,{TAU(J),J=1,JJ)

34 6 FORMATC МОМЕШи ВРЕМЕНИ :/(5G11.3))

35 С ЗАПИСЬ СОБСТВЕННЫХ ЧИСЕЛ В МАССИВ А

36 BIO=Alf«DL/ALMl

37 CALL KORH1(HN,0.01,BIO.A)

36 С ВЫЧИСЛЕНИЕ ЭЛЕМЕНТОВ РАБОЧИХ МАССИВОВ В1 Я В?

39 Ю 7 N=.l,NN

40 AN=A(N)

41 С КВАДРАТ ИОРШ СОБСТВЕННОЙ ФУШЩ

42 F2=DL*(l.+BI0/(BI0i(BIO+AN4AN))/2.

43 С КОМПЛЕКС B1(N)

44 B1(N)=T0«DL»SIN(AH)/(AN*F2)

45 С ИЗОБРАШЕНИЕ РАСПРЕДЕЛЕНИЯ МОЩНОСТИ ВНУТР. ИГГСЧНИКОВ Ш

46 CALL SIMPS(0.,DL,M,WINT.llW)

47 с КОМПЛЕКС B2(N)

48 7 B2(N)-WN»DL*DL/(ALM1#AN»F2)

49 С ЦИКЛ ПО МОМЕНТАМ ВРЕМЕНИ

50 DO И J-l,jJ

51 С ВЫЧИСЛЕНИЕ КОМПЛЕКСОВ С(Ы>

52 DO S N-1,NN

53 AN-A{N)

Рис. 2.12

55 56 57 56 59 60 81 62 63 64 65

86 70 71

72 73 74 75 78 77

79 79 60 91 82 63 64 65 86 67

90 91 92 93 94 95 96 97 88 99 190 101 102 103 104 105 106

Зак. 2217

FO=AT»TAU(J)/(DLi.DL)

8 C(N)=EXP{-ANxAN«FO)

1ЩКЛ ПО ПРОСТРАНСТВЕННЫМ ТОЧКАМ DO 9 1=1,11

СУММИРОВАНИЕ РЯДА, ВЫЧИСЛЕНИЕ ПЕРЕГРЕВОВ Т(1)

Т(1)=0.

DO В N=1.NN

9 T(I)=T{I)+{BI(N)»C(N)+BZ(N)»(1.-C(N)))*C0S(A(W1#X<I)/D.) ПЕЧАТЬ ТЕМПЕРАТУР Б ДАННЫЙ МОМЕНТ ВРЕМЕНИ

PRINT 10,TAU(J),(T(I),1=1.111 10 FORMAT(/ ВРЕМЯ -,G11.3/ ТЕМПЕРАТУРЫ :7(6011-3)) И CONTINUE

STOP

ПОДПРОГРАММА ПОДЫНТЕГРАЛЬНОЙ ФУНКЦИИ

ДЛЯ ВЫЧИСЛЕНИЯ ИЗОБРАЖЕНИЯ РАСПРВДЕЛЕНйЯ МОЩНОСТМ

FUNCTION WINTCX)

COMMON /INT/ AN,DL

WINT=0V(X)«COS(AN»X/DL)

RETURN

ПОДПРОГРАММА-ФУНКЦИЯ РАСПРЕДЕЛЕНИЯ МОЩНОСТО

ВНУТРЕННИХ ИСТОЧНИКОВ

FUNCTION QVCX)

QV0=1.E6

DL=0.01

0V=QV0» EXP(-(3»X/DL)»»г)

RETURN

ПОДПРОГРАММА ВЫЧИСЛЕНИЯ КОРНЕЙ ТРАНСЦИЩЕНШОГО УРАВНЕНИЯ CTG(X)=X/BIO МЕТОДОМ ПОЛОВИННОГО ДЕЛЕНИЯ SUBROUTINE KORNI «(N,EPS,B10,X) ВХОДНЫЕ ПАРАМЕТРЫ N-ЧИСЛО КОРНЕЙ

ЕРЗ-ПОГРЕШНОСГЬ ОПРЕДЕЛЕНИЯ КОРНЯ ВЮ-ЧИСЛО БИО

ВЫХОДНОЙ ПАРАНЕ!? : X(N)-HACCHB КОРНЕЙ

DIMENSION Х{1) F(X)=X/B10-C0S(X)/SIN(X) НАЧАЛЬНЫЙ ИНТетВАЛ ДЛЯ ПЕРВОГО КОРНЯ 1=1

В=1.57Й79

ДЕЛЕНИЕ ИНТЕРВАЛА ПОПОЛАМ 1 С=(А+В)/2.

ВЫЧИСЛЕНИЕ ФУНКЦИИ В ТОЧКЕ ДЕЛЕНИЯ FC=F(C)

ТЧЛБОР НОВОЙ ГРАНИЦЫ ИНТЕРВАЛА

Рис. 2.12. Продолжение