11. Приближенное вычисление интеграла по формуле трапеций

Постановка задачи.

Численное интегрирование функции сводится к непосредственному расчету определенного интеграла по ряду значений подынтегральной функции. Типичная методика заключается в том, что заданную функцию f(x) на сегменте [а, b] заменяют интерполирующей или аппроксимирующей функцией f(х) простого вида, полагая

(2.37)

Заметим, что функция f(х) выбирается такой, чтобы соответствующий интеграл вычислялся непосредственно.

Формула трапеций.

Найдем коэффициенты Котеса при п = 1, т.е.

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

2_2

т.е.

(2.44)

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

(2.45)

На практике используют более удобную формулу расчета погрешности

(2.46)

где x - некоторое число в интервале (х₀, х₁).

Заметим, что при у" > 0 по формуле (2.44) вычисляют интеграл с избытком, а при у" < 0 - с недостатком.

12. Приближенное вычисление интеграла по формуле Симпсона

2_2

Формула Симпсона.

Вновь обращаясь к (2.42), найдем коэффициенты Котеса при п = 1. Тогда имеем

Н₀ = 1/6, Н₁= 2/3, Н₂ = 1/6.

Поскольку x₂-x₀ = 2h, то для рассматриваемого случая трех точек (рис. 2.2)

из (2.43) получаем формулу Симпсона

(2.47)

Как видно из рис. 2.2, формула Симпсона получается в результате замены данной кривой у =f(х) параболой у = L₂(x) (показана штриховой линией), проходящей через три точки A₀, A₁ и A₂.

В этом случае для погрешности расчета находят

(2.48)

Откуда после несложных преобразований получают в более удобном виде:

(2.49)

где x - некоторое число в интервале [х₀, х₂].

Анализ формулы (2.49) свидетельствует о том, что формула Симпсона является точной для многочленов как второй, так и третьей "степени, поскольку в этом случае y¹^V(x) = 0 и имеем нулевую погрешность при расчете (R = 0).

На практике при машинном расчете определенных интегралов удобно погрешность вычислений оценивать по следующей методике. Сначала вычисляют значение интеграла I_h, при заданном шаге, а затем находят новое значение I₂_h, при вдвое большем шаге. После этого находят абсолютную погрешность

13. Матричная алгебра. Методы расчета определителя матрицы

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

Каждый элемент матрицы снабжают двумя индексами, первый из которых указывает номер его строки, а второй - номер столбца. Например, элемент a_lk располагается в l-й строке и k-м столбце соответствующей матрицы. Матрицу A, составленную из a_ij элементов, можно записать различными способами, например

(5.1)

Матрицы характеризуются размером, под которым понимают число его строк т и число столбцов п, т.е. т x п. С этой точки зрения различают прямоугольные (т¹ п) и квадратные (т=п) матрицы. К прямоугольным матрицам относят также матрицы-столбцы (п = 1) или матрицы-строки (т = 1). Эти матрицы называют также векторами, приписывая каждому из элементов такой матрицы роль соответствующей компоненты вектора. Квадратную матрицу размера п x п относят к матрице п -го порядка.

Различают следующие разновидности квадратных матриц ||a_ij || порядка п:

· верхняя треугольная матрица, если a_ij= 0 для всех i >j;

· нижняя треугольная матрица, если a_ij= 0 для всех i <j;

· диагональная матрица, если a_ij= 0, для всех i ¹ j;

· единичная матрица (частный случай диагональной), для которой имеют a_ij= d_ij, где d_ij - символ Кронекера, т.е. d_ij = 0 при i ¹ j и d_ij = 1 при i= j.

Для конкретизации указанных разновидностей ниже представлены соответствующие матрицы, например, третьего порядка:

Матрица, все элементы которой равны нулю, называется нулевой.

В квадратных матрицах различают две диагонали: главную (проходит из левого верхнего угла в правый нижний) и

побочную (проходит из правого верхнего угла в левый нижний).

Матрица, симметричная относительно главной диагонали, называется симметрической.

Переставляя в квадратной матрице

(5.3)

строки со столбцами, получают транспонированную матрицу, т.е.

(5.4)

Определителем (детерминантом) квадратной матрицы n-го порядка называют число D, образованное из n² ее элементов (чисел) следующим образом:

где a,b,…q — всевозможные перестановки для набора чисел из п сомножителей в каждом, содержащие по одному элементу из каждой строки и столбца;

k — коэффициент учета знака (k = 0 при четном числе перестановок и k = 1 при нечетном).

Определитель можно обозначить также символами

(5.6)

и т.д.

В общем случае отсутствия нулевых элементов в матрице n-го порядка ее определитель (5.5) будет содержать n! (n-факториал) слагаемых, каждый из которых имеет п сомножителей .

Матрица, определитель которой равен нулю, называется особенной или вырожденной. Соответственно при det||A||¹0 получают неособенную (невырожденную) матрицу.

Минором М_ij элемента a_ij в матрице n-го порядка называют определитель (п - 1)-го порядка, получаемый из исходной матрицы вычеркиванием (i-й строки j-го столбца.

Соответственно алгебраическим дополнением А_ij элемента a_ij,у именуют произведение минора М_ij на

(-1)^i+j:

A_ij=(-1)^i+j М_ij (5.7)

Рангом матрицы ||A|| называют целое число р, равное наибольшему порядку ее неособенной квадратной подматрицы. В случае rang||A||=р любые другие подматрицы более высокого порядка (выше р) должны быть вырожденными (особенными).

Непосредственный расчет определителей по формуле (5.5) не имеет практического смысла для матриц порядка n>3 из-за резкого возрастания числа слагаемых. Поэтому на практике используют различные искусственные приемы снижения трудоемкости вычислений. Они основаны на простейших преобразованиях исходной матрицы к более удобному для расчета определителя виду. К ним можно отнести, например, верхнюю треугольную матрицу (5.2), определитель которой равен произведению ее диагональных элементов.

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

а) при перестановке двух любых строк (столбцов) определитель умножается на -1;

б) общий множитель элементов любой строки (столбца) можно выносить за знак определителя;

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

Процедуру приведения исходной матрицы к треугольному виду с регулярным использованием указанных выше преобразований относят к методу исключения Гаусса.

14. Методы расчета обратной матрицы. Обращение матриц методом преобразование единичной матрицы

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

Обозначив матрицу, обратную||A||, через ||A||^-1, по определению получим

||A||•||A||^-1=||A||^-1•||A||=||E||, (5.31)

где ||E|| - единичная матрица.

Процедуру нахождения обратной матрицы

||A||^-1называют обращением матрицы.

Прямой способ обращения матриц.

В этом случае обратную матрицу вычисляют по формуле

_(5.32)

где |A| — определитель исходной матрицы; А_ij— алгебраическое дополнение элемента a_ij (5.7) в матрице ||A||.

Как видно из (5.32), обратные матрицы существуют только для неособенных матриц, определитель которых не равен нулю. Обратная матрица, если она существует, является единственной.

Метод получения обратной матрицы через алгебраические дополнения элементов исходной матрицы (5.32) не находит применения в ЭВМ из-за чрезвычайно большого числа вычислений для матриц высоких порядков.

Косвенный способ обращения матриц.

Этот способ, широко используемый при расчете обратной матрицы в ЭВМ, заключается в следующем. Приписывают к исходной матрице ||A|| справа единичную матрицу ||E||. Выполняя известные преобразования (см. разд. 5.2) над cоставной матрицей, приводят подматрицу ||E|| к единичному типу ||A||^-1. Иногда на месте первоначально приписанной подматрицы ||E|| получают искомую обратную матрицу ||A||^-1, т.е.

(5.33)

или соответственно

_(5.34)

где матрица ||E|| квадратная п-го порядка, а стрелки обозначают соответствие между исходной составной матрицей

и преобразованной

;

штриховая вертикальная линия в составных

матрицах выполняет роль разделительной черты.

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

Каждый столбец подматрицы ||A|| в составной матрице

подвергают преобразованию (обработке), начиная с первого и кончая n-м. Процедуру такой обработки столбцов, рассматриваемую в общем виде, можно алгоритмизировать следующим образом.

1. Обозначим очередной столбец подматрицы ||A||, подлежащий обработке, через j-й.

Пусть в j-м столбце

H_ij=|a_ij|_max, j£ i£ n, (5.35)

Является элементом наибольшим по абсолютной величине.

Тогда выполняют следующее:

а) если H_ij = 0, то подматрица вырожденная и ее обращение невозможно (процедура на этом заканчивается);

б) если H_ij¹ 0, то нужно в составной матрице

поменять местами i-ю строку (с H_ij) на j-ю, обеспечивая размещение элемента с наибольшей абсолютной величиной на месте диагонального элемента a_ij .

2. Каждый элемент j-й строки матрицы

делят на значение a_ij = H_ij получая при этом на месте диагонального элемента единицу. После этого для всех i-х строк (i¹j) матрицы

выполняют поэлементные вычитания вида:

в подматрице ||A||

(5.36)

где k - номер элемента в i-й строке (1£ k £ п), а знак "штрих" служит для обозначения полученного (нового) значения элементов в результате операции (5.36);

в подматрице ||A||

(5.37)

где е - элемент единичной подматрицы ||E||.

3. В результате выполнения указанных операций составная матрица

преобразуется к эквивалентной матрице

из которой легко записать значения всех элементов искомой обратной матрицы [A]^-1, как это показано в (5.34).

15. Собственные значения матрицы. Обращение матриц методом преобразование единичной матрицы

Собственные векторы матрицы.

Пусть ||A|| - квадратная матрица п-го порядка, а l - некоторое число.

Любой ненулевой столбцовый вектор ||X||, для которого

||A|| ||X||= l||X||, (5.38)

называют собственным вектором матрицы ||A||, а число l - собственным значением(характеристическим числом) этой матрицы.

Характеристическое уравнение матрицы.

Равенство (5.38) эквивалентно уравнению

(||A||- l||E||) ||X||=0 (5.39)

которое отображает в матричной форме систему п линейных однородных уравнений относительно компонентов собственного вектора ||X||, т.е.

(a₁₁ - l)x₁ + a₁₂x₂+ … + a_1n x_n = 0,

a₂₁x₁ + (a₂₂ - l)x₂+ … + a_2n x_n = 0, (5.40)

………………………………………

a_n1x₁ + a_n2x₂+ … + (a_nn - l) x_n = 0.

Уравнения (5.40) имеет ненулевые решения тогда и только тогда, когда их определитель равен нулю, т.е. при условии

_(5.41)

или сокращенно

det(||A||- l||E||) = 0. (5.42)

Уравнения (5.41) или (5.42) можно представить и в виде многочлена

lⁿ - s₁lⁿ^-¹ + s₂lⁿ^-²+ ... + (-l)ⁿ^-1s_n_-1l + (-1)ⁿs_n=0. (5.43)

Многочлен, стоящий в левой части (5.43), называют характеристическим, а само уравнение - характеристическим уравнением матрицы ||A||. Таким образом, уравнения (5.41)-(5.43) являются разновидностями одного и того же характеристического уравнения, решение которого обеспечит определение собственных значений l_k, k=1,2,...,п, матрицы ||A|| порядка п.

Заметим, что каждому собственному значению l_k матрицы ||A|| соответствует совокупность из п значений х_i, i=1,2,...,п, являющихся компонентами ее собственного вектора ||X|| для заданного l_k. Эти компоненты находят из решения (5.40) при значении

l=l_k.

Собственные значения симметрических матриц.

При расчете реальных устройств, как правило, получают симметрические матрицы с вещественными элементами. Такие системы характеризуются следующими специфическими свойствами, упрощающими анализ и решение характеристического уравнения (5.42):

а) если матрица ||A|| симметрическая, то для любых двух ее векторов ||X|| и ||Y|| с компонентами х₁, х₂,..., х_n и y₁, y₂,..., y_n соответственно можно записать скалярные произведения

(||X||,||A||• ||Y||) = (||A||•||X||,||Y||); (5.44)

б) собственные значения вещественной

симметрической матрицы вещественны;

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

Методы расчета собственных значений.

Вопрос определения собственных векторов ||X|| исходной матрицы ||A|| чрезвычайно сложен. Как отмечалось, для этого требуется реализация следующих процедур:

а) от исходной матрицы ||A||, например, п-го порядка перейти к уравнению (5.42) и на этой основе получить характеристический многочлен типа (5.43);

б) решить характеристическое уравнение (5.43) п-го порядка и получить численно все собственные значения l_i, i = 1, 2,..., п, исходной матрицы ||A||;

в) подставив каждое из i-x значений l_i в (5.40), решить эту систему уравнений для определения соответствующих компонентов собственных векторов ||X||_i, i = 1, 2,..., п.

Для исходных матриц ||A|| высокого порядка непосредственная реализация указанных процедур требует громоздких вычислений. В то же время собственные значения l₁, l₂,..., l_n матрицы ||A|| позволяют получить решение дифференциальных уравнений соответствующего порядка. Поэтому разработаны специальные косвенные способы решения задачи, заметно упрощающие трудоемкость вычислений.

Применительно к симметрическим исходным матрицам ||A|| особую актуальность для реализации на ЭВМ имеет итерационный метод численного расчета собственных значений, разработанный Якоби. Он основан на применении специальных преобразований подобия с использованием ортогональной матрицы элементарных вращений Якоби, обозначаемой ||Т||.

Для понимания преобразований подобия в методе Якоби ниже приводится матрица элементарных вращений:

Процесс преобразований подобия состоит в построении последовательности матриц

||A||=||A'||,||A''||,||A'''||,..., каждая из которых получается из предыдущей при помощи элементарного шага (вращения). Заметим, что матрица вращения (5.45) выбирается на каждом шаге так, чтобы элемент р-й строки и q-го столбца преобразованной матрицы обращался в нуль. При этом если пары индексов (р, q) выбираются так, что на каждом шаге аннулируется наибольший по абсолютной величине внедиагональный элемент преобразованной матрицы, то вычислительный процесс сходится.

Если исходная симметрическая матрица ||A|| не имеет кратных собственных значений, то последовательное аннулирование всех ее внедиагональных элементов (с точностью до e) приведет к тому, что ее диагональные элементы приблизятся к собственным значениям l_i, (с точностью до e²).

Реализация метода Якоби начинается с выявления наибольшего по абсолютной величине внедиагонального элемента в исходной симметрической матрице ||A||, например:

а_рq = а_qр. (5.46)

Соответствующий угол поворота j, обращающий эти элементы в нуль, находят по формуле

ctg2j = (а_рр – a_qq)/2a_pq, -p/4<j<p/4, (5.47)

где а_рр и a_qq - соответствующие диагональные элементы исходной матрицы ||A||.

Тогда в преобразованной матрице

||A'||=||Т||^Т_pq•||A||•||Т||_pq, (5.48)

на месте элементов а_р_q и a_qp, будут стоять нули.

Вновь применяя к матрице ||A'|| элементарное вращение, можно в матрице ||A''|| обратить в нуль еще два новых элемента. При этом нулевые элементы матрицы |A'| остаются без изменений и т.д.

В общем виде (для k-го шага) алгоритм рассмотренного метода Якоби следующий:

а) выбор элемента а⁽^k⁾_р_q из матрицы ||A⁽^k⁾||;

б) вычисление по итерационной схеме

||A⁽^k⁺¹⁾||=||Т||^Т_pq•||A⁽^k⁾||•||Т||_pq, k=1,2,...; (5.49)

в) окончание (обрыв) расчета, если

e>0 и задано, (5.50)

в противном случае - переход к (k + 1)-му шагу.

При обрыве матрица ||A⁽^k⁾|| имеет диагональный вид. С точностью, определяемой e, диагональные элементы являются приближениями собственных значений l_i, а столбцы матрицы произведения всех применяемых ||Т||_pq — приближениями собственных векторов.

Сходимость итерационного процесса обеспечивается при выборе в качестве ведущего элемента а⁽^k⁾_р_q наибольшего по абсолютной величине внедиагонального элемента, т.е. при условии

а⁽^k⁾_р_q =| а⁽^k⁾_р_q |_max ,

16. Численное решение систем линейных уравнений

Классы задач линейной алгебры

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

Ø решение систем линейных алгебраических уравнений (СЛАУ);

Ø вычисление определителей матриц;

Ø нахождение обратных матриц;

Ø определение собственных значений и собственных векторов матриц

Постановка задачи решения СЛАУ: (1), где – квадратная матрица коэффициентов размерности n, – вектор неизвестных, – вектор свободных коэффициентов. Иногда СЛАУ представляют в виде расширенной матрицы размерности n×(n+1), где в качестве последнего столбца фигурирует вектор свободных коэффициентов. В координатном представлении такая СЛАУ выглядит следующим образом:

Для решения СЛАУ применяют в основном два класса методов: прямые (выполняемые за заранее известное количество действий) и итерационные (обеспечивающие постепенную сходимость к корню уравнения, зависящую от многих факторов). Прямые методы обычно применяются для решения систем порядка n < 200, для бóльших n используются итерационные методы. Перед решением СЛАУ требуется проанализировать корректную постановку задачи:

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

Если не имеет элементов с большими по модулю значениями – решение устойчиво (см. пример к главе 1). Показателем плохо обусловленных систем является

19. Простейшие численные методы интегрирования дифференциальных уравнений (Методы Эйлера)

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

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

Вывод формулы.

Пусть дано дифференциальное уравнение

(3.21)

и начальное условие

y(х0)=у0. (3.22)

1 Требуется найти решение уравнения (3.21) методом Эйлера, т.е. получить интегральную кривую у(х) на отрезке [а, b]. Выбрав достаточно малый шаг h, находим число интервалов на отрезке [а, b], т.е.

п = (b - a)/h. (8.23)

3_2

Разделив отрезок [а, b] на и одинаковых интервалов шириной h (см. рис. 3.2),

Рис 3.2

получаем систему равноотстоящих точек

хi = х0 + ih, i = 0,1,2,.... п, (3.24)

где x0 =а; хn= b.

Полагая на i-м интервале функцию f(хi,уi)=у'(хi) прямолинейной, а ее крутизну tgai, = у'(xi), исходное уравнение (3.21) можно представить (рис. 3.3) как

(3.25)

Обозначив хi+1- хi =h и f(хi,уi)=у'(хi), получим итерационную формулу •Эйлера, т.е.

yi+1=yi+hy'i. (3.26)

Применяя рассмотренный метод Эйлера сначала для первого интервала на отрезке [а, b], а затем для второго и т.д., можно получить последовательно значения ординат уi, i = 1,2,..., п, определяющие интегральную кривую решения (8.21) в виде ломаной Эйлера, показанной на рис. 3.4.

Анализ результатов.

Рассмотренный явный метод Эйлера (3.26) обладает низкой точностью и способностью накапливать погрешность.

Сопоставляя (3.26) с (3.10), можно заключить, что метод Эйлера является частным случаем решения посредством ряда Тейлора, в котором используют лишь первые два члена ряда, т.е. имеем порядок п = 1.

Как видно из (3.12), формула Эйлера может обеспечить на каждом шаге точность до величин порядка h2. Поэтому для суммарной точности метода на отрезке [а, b] имеем величину nh2, пропорциональную шагу h, поскольку значение п обратно пропорционально h (3.23). В связи с этим при расчетах требуется выбирать значение шага h, обеспечивающее заданную суммарную точность e или точность на шаге ei.

Пусть, например, требуется получить точность расчета на шаге ei =0,01. Тогда значение шага будет:

Обращаясь кc(3.23), находим число интервалов п, требуемое для обеспечения заданной точности.

Приближенную оценку погрешности на отрезке [а, b] при выбранном шаге h можно получить также по значению контрольной разности

где уh и у2h - результаты вычислений, полученные соответственно при шаге h и 2h.

Неявный метод Эйлера основан на следующем предположении. Если в каждом i-м шаге расчетной формулы типа (3.26) использовать значение производной не в начальной точке, а в конечной, то имеют

yi+1=yi+hy'i. (3.28)

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

20. Методы Рунге-Кутты

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

Обоснование методов Рунге-Кутта

Для исходного дифференциального уравнения у'=f(x,у) определяют на отрезке x₀£x£x₀+h точки (x_i,h_i) и константы R_i

i = 0, 1, 2,..., п. При этом ставится задача, чтобы интегральная кривая решения у =f₁(x), проходя через точку (x₀,y₀) имела ординату у₁=f₁(x₀+h), возможно точнее совпадающую с выражением

k=y₀+h[R₀ f(x₀,h₀)+...+ R_i f(x_i,h_i)+..+ R_n f(x_n,h_n)]. (3.34)

Термин "возможно точнее" означает, что при разложении аппроксимирующей f(x) и исходной f(x,y) функций в ряд по степеням п в правой и левой частях равенства

f₁(x₀+h) = f[xo + h, f(x₀ + h)]

должно совпадать наибольшее число членов.

При этом должны быть обеспечены следующие соотношения:

x₀=x₀	h₀=y₀
x₁= x₀+a₀h	h₁= y₀+b₀k₁	k₁= h f(x₀,h₀);
x₂= x₀+a₁h	h₂= y₀+b₁k₁+g₁k₂	k₂= h f(x₁,h₁);

Заметим, что R_i,a_i,b_i…. выбираются независимо от функции f(x,y).

Методика решения.

Для интегральной кривой y=f(x), проходящей через точку (х₀, у₀) в равноотстоящих друг от друга точках x₀, x₁, x₂,..(x_i₊₁-x_i=h>0) получают такие значения y₀, y₁, y₂... чтобы обеспечить у_i»f_i(x) .При этом в зависимости от заданной точности для последовательного определения y₁, y₂.... используют формулы Рунге-Кутта соответствующих порядков, приведенные в табл. 3.1 и 3.2.

Таблица 3.1 1 вариант формул Ронге-Кутта различных порядков точности

Порядок точности	Расчетное соотношение	Значение коэффициентов
Второй	y_i+1= y_i+(k₁+ k₂)/2	k₁=hf(x_i,y_i) k₂=hf(x_i+h,y_i+k₁)
Третий	y_i+1= y_i+k₁/6+ 2k₂/3+ k₃/6	k₁=hf(x_i,y_i) k2=hf(x_i+h/2,y_i+k₁/2) k₃=hf(x_i+h,y_i+k₁+2k₂)
Четвертый	y_i₊₁₌y_i+k₁/6+k₂/3+k₃/3+k₄/6	k₁=hf(x_i,y_i) k₂=hf(x_i+h/2,y_i+k₁/2) k₃=hf(x_i+h/2,y_i+k₂/2) k₄=hf(x_i+h, y_i+k₃)

Таблица 3.2 2 вариант формул Ронге-Кутта различных порядков точности

Порядок точности	Расчетное соотношение	Значение коэффициентов
Второй	y_i+1= y_i+k₂	k₁=hf(x_i,y_i) k₂=hf(x_i+h/2,y_i+k₁/2)
Третий	y_i₊₁₌ y_i+k₁/4+ 3k₃/4	k₁=hf(x_i,y_i) k2=hf(x_i+h/3,y_i+k₁/3) k₃=hf(x_i+2h/3,y_i+k₁+2k₂/3)
Четвертый	y_i₊₁₌y_i+k₁/8+3k₂/8+3k₃/8+k₄/8	k₁=hf(x_i,y_i) k₂=hf(x_i+h/3,y_i+k₁/3) k₃=hf(x_i+2h/3,y_i+2k₂/3) k₄=hf(x_i+h, y_i+k₁-k₂+k₃)

Практические рекомендации.

Для приблизительной оценки погрешности при расчетах на ЭВМ удобно проводить вычисления сначала с шагом h, а затем - 0,5h, получая разность между этими результатами e_y

После этого находят погрешность приближенного решения для формулы n-го порядка, пользуясь соотношением

e_p=|e_y|/(2ⁿ-1) где e_y=y_h-y_h_/2

21. Метод прогонки

Прогонкой называется модификация метода Гаусса для решения систем линейных алгебраических уравнений с трехдиагональной матрицей. Если матрица системы обладает определенными свойствами, то метод прогонки является численно устойчивым и очень эффективным методом, который позволяет практически мгновенно решать одномерные краевые задачи, одну из которых мы рассмотрели в предыдущем разделе. Большинство корректно поставленных физических задач приводит к системе уравнений с хорошей матрицей, и в этих случаях метод прогонки проявляет слабую чувствительность как к погрешностям задания начальных условий, так и к погрешностям вычислительного характера.

В предыдущем разделе была сформулирована так называемая первая краевая задача, в которой требуется найти значения функции во внутренних узлах сетки при условии, что на границах области они известны. В теории и на практике рассматриваются задачи с более сложными граничными условиями. Например, когда на одной из границ известна не функция, а ее первая производная — граничное условие второго рода. Имеют место и постановки задач с граничными условиями третьего рода, когда на границе должно выполняться какое-то известное заранее соотношение между функцией и ее первой производной. С точки зрения численной реализации все три типа задач можно описать с помощью соотношений одного и того же вида:

U0=y0U1+б0, (6)

Un=ynUn-1+бn, (7)

Они связывают значения разностных аналогов Ui, непрерывной функции U(x) в двух узлах, прилегающих к левой или правой границе. Так, граничное условие первого рода иUo = с может быть задано с помощью пары параметров: у0= 0, б0 = с, а условие второго рода dU/dx|0= с с помощью другой лары: у0 = 1,бo=ch, где h — это шаг сетки. В нашем приложении будет работать немодальный диалог, который позволит пользователю задавать различные типы граничных условий, изменяя численные значения четырех коэффициентов уo, бo, yn, бn

Суть метода прогонки заключается в том, что, используя специфику структуры матрицы системы уравнений (наличие трех диагоналей), удается получить рекуррентные формулы для вычисления последовательности коэффициентов прогонки, которые позволяют на обратном ходу вычислить значения функции в узлах сетки. Рассматривая конечно-разностное уравнение для первой тройки узлов:

b1U1+c1U2=-a1U0,

видим, что оно совпадает по форме с обобщенным граничным условием (6) и связывает между собой два соседних значения U1, и U2. Перепишем его в виде:

d1U2+e=U1, (8)

где d1 и е1вычисляются по известным значениям. Наблюдательный читатель заметит, что это справедливо только для задач первого рода. Чуть позже мы получим общее решение. Теперь мы можем исключить £/, из уравнения для следующей тройки узлов:

a2U1+b2U2+c2U2=f2,

подставив значение U1 из уравнения (8). После этой процедуры последнее уравнение также может быть приведено к виду:

d3U3+e2=U2,

Подстановки можно продолжать и дальше, но для получения рекуррентного соотношения, достаточно рассмотреть одну из них для произвольного индекса i. Подставив

di-1Ui+ei-1=Ui-1,

в уравнение

aiUi-1+biUi+ciUi+1=fi,

получим:

Ui=-[CiUi+1/(aidi-1+bi)]+[fi-ai+1*ei+1/(aidi-1+bi)] (9)

Это соотношение дает две рекуррентные формулы для коэффициентов:

di=-Ci/(ai*di-1+bi) (10)

ei=(fi-ai*ei-1)/(aidi-1+bi) (11)

Цикл вычисления последовательности коэффициентов в соответствии с этими формулами носит название прямого хода прогонки. Начальные значения коэффициентов определяются предварительно заданным граничным условием (6):

do=yo, eo=бo,

Цикл прямого хода повторяется N-1 раз. Последними будут вычислены коэффициенты dN-1 и eN-1, которые связывают функции в двух узлах вблизи правой границы:

Un-1=dn-1Un+en-1 (12)

Если на правой границе задано условие первого рода Un = с, то уже можно вычислить Un-1 по формуле (12) и далее продолжать обратный ход прогонки, используя уравнения (9) при I = N - 1,..., 1, 0. Если условие более сложное, то вместе с уравнением (12) надо рассмотреть уравнение (7), определяющее граничное условие на правой границе. Напомним его:

Un=ynUn-1+бn (7)

Соотношения (7) и (12) составляют систему из двух уравнений с двумя неизвестными. Используя определители, запишем ее решение.

Un-1=(en-1+бndn-1)/(1-yndn-1) (13)

Un=(бn+ynen-1)/(1-yndn-1)

Таким образом, мы нашли значения в двух узлах, лежащих вблизи правой границы расчетной области. Теперь, используя формулу (9) и уменьшая индекс i от N= 2 до 0, можно вычислить все неизвестные £/.. Этот процесс носит название обратного хода прогонки. Почему-то в голову приходит лозунг нашего времени: «Цели ясны, задачи определены. За работу, товарищи!» Нам осталось всего лишь реализовать описанный алгоритм в виде MFC-приложения.

17. Метод Гаусса.

Ме́тод Га́усса— классический метод решения системы линейных алгебраических уравнений (СЛАУ). Это метод последовательного исключения переменных, когда с помощью элементарных преобразований система уравнений приводится к равносильной системе ступенчатого (или треугольного) вида, из которого последовательно, начиная с последних (по номеру) переменных, находятся все остальные переменные.

Пусть исходная система выглядит следующим образом

$\left\{\begin{array}{lcr} a_{11}x_1+\ldots+a_{1n}x_n &=& b_1 \\ \ldots & & \\ a_{m1}x_1+\ldots+a_{mn}x_n &=& b_m \\ \end{array}\right.\iff Ax = b,\quad A=\left( \begin{array}{ccc} a_{11} & \ldots & a_{1n}\\ \ldots & & \\ a_{m1} & \ldots & a_{mn} \end{array}\right),\quad x = \left( \begin{array}{c}x_1 \\ \vdots \\ x_n \end{array} \right), \quad b = \left( \begin{array}{c}b_1 \\ \vdots \\ b_m \end{array} \right).\quad (1)$

Матрица $A$ называется основной матрицей системы, $b$ — столбцом свободных членов.

Тогда согласно свойству элементарных преобразований над строками основную матрицу этой системы можно привести к ступенчатому виду(эти же преобразования нужно применять к столбцу свободных членов):

$\left\{\begin{array}{rcl} \alpha_{1j_1}x_{j_1}+\alpha_{1j_2}x_{j_2}+\ldots+\alpha_{1j_r}x_{j_r}+\ldots+\alpha_{1j_n}x_{j_n} &=& \beta_1 \\ \alpha_{2j_2}x_{j_2}+\ldots+\alpha_{2j_r}x_{j_r}+\ldots+\alpha_{2j_n}x_{j_n} &=& \beta_2 \\ &\ldots& \\ \alpha_{rj_r}x_{j_r}+\ldots+\alpha_{rj_n}x_{j_n} &=& \beta_r \\ 0 &=& \beta_{r+1} \\ &\ldots& \\ 0 &=& \beta_m \end{array}\right.,\qquad \alpha_{1j_1},\ldots,\alpha_{rj_r}\neq 0.$

При этом будем считать, что базисный минор (ненулевой минор максимального порядка) основной матрицы находится в верхнем левом углу, то есть в него входят только коэффициенты при переменных $x_{j_1},\ldots,x_{j_r}\!$ .

Тогда переменные $x_{j_1},\ldots,x_{j_r}\!$ называются главными переменными. Все остальные называются свободными.

Если хотя бы одно число $\beta_i\neq 0\!$ , где $i > r$ , то рассматриваемая система несовместна.

Пусть $\beta_i = 0\!$ для любых $i > r$ .

Перенесём свободные переменные за знаки равенств и поделим каждое из уравнений системы на свой коэффициент при самом левом $x\!$ ( $\alpha_{ij_i},\, i=1,\ldots,r\!$ , где $i\!$ — номер строки):

$\left\{\begin{array}{rcc} x_{j_1}+\widehat{\alpha}_{1j_2}x_{j_2}+\ldots+\widehat{\alpha}_{1j_r}x_{j_r}&=& \widehat{\beta}_1-\widehat{\alpha}_{1j_{r+1}}x_{j_{r+1}}-\ldots- \widehat{\alpha}_{1j_n}x_{j_n} \\ x_{j_2}+\ldots+\widehat{\alpha}_{2j_r}x_{j_r}&=& \widehat{\beta}_2-\widehat{\alpha}_{2j_{r+1}}x_{j_{r+1}}-\ldots- \widehat{\alpha}_{2j_n}x_{j_n} \\ &\ldots& \\ x_{j_r}&=& \widehat{\beta}_r-\widehat{\alpha}_{rj_{r+1}}x_{j_{r+1}}-\ldots- \widehat{\alpha}_{rj_n}x_{j_n} \\ \end{array}\right., \qquad \widehat{\beta}_i=\frac{\beta_i}{\alpha_{ij_i}},\quad \widehat{\alpha}_{ij_k}=\frac{\alpha_{ij_k}}{\alpha_{ij_i}}\quad (2)$ ,
где $i=1,\ldots,r,\quad k=i+1,\ldots,n.\!$

Если свободным переменным системы (2) придавать все возможные значения и решать новую систему относительно главных неизвестных снизу вверх (то есть от нижнего уравнения к верхнему), то мы получим все решения этой СЛАУ. Так как эта система получена путём элементарных преобразований над исходной системой (1), то по теореме об эквивалентности при элементарных преобразованиях системы (1) и (2) эквивалентны, то есть множества их решений совпадают.

Алгоритм

Алгоритм решения СЛАУ методом Гаусса подразделяется на два этапа.

На первом этапе осуществляется так называемый прямой ход, когда путём элементарных преобразований над строками систему приводят к ступенчатой или треугольной форме, либо устанавливают, что система несовместна. А именно, среди элементов первого столбца матрицы выбирают ненулевой, перемещают его на крайнее верхнее положение перестановкой строк и вычитают получившуюся после перестановки первую строку из остальных строк, домножив её на величину, равную отношению первого элемента каждой из этих строк к первому элементу первой строки, обнуляя тем самым столбец под ним. После того, как указанные преобразования были совершены, первую строку и первый столбец мысленно вычёркивают и продолжают пока не останется матрица нулевого размера. Если на какой-то из итераций среди элементов первого столбца не нашёлся ненулевой, то переходят к следующему столбцу и проделывают аналогичную операцию.
На втором этапе осуществляется так называемый обратный ход, суть которого заключается в том, чтобы выразить все получившиеся базисные переменные через небазисные и построить фундаментальную систему решений, либо, если все переменные являются базисными, то выразить в численном виде единственное решение системы линейных уравнений. Эта процедура начинается с последнего уравнения, из которого выражают соответствующую базисную переменную (а она там всего одна) и подставляют в предыдущие уравнения, и так далее, поднимаясь по «ступенькам» наверх. Каждой строчке соответствует ровно одна базисная переменная, поэтому на каждом шаге, кроме последнего (самого верхнего), ситуация в точности повторяет случай последней строки.

Метод Гаусса требует порядка $O (n 3)$ действий.

Простейший случай

В простейшем случае алгоритм выглядит так:
$\left\{\begin{array}{lcc} a_{11} \cdot x_1 + a_{12} \cdot x_2 + \ldots + a_{1n} \cdot x_n & = b_1 & (1) \\ a_{21} \cdot x_1 + a_{22} \cdot x_2 + \ldots + a_{2n} \cdot x_n & = b_2 & (2) \\ \ldots & & \\ a_{m1} \cdot x_1 + a_{m2} \cdot x_2 + \ldots + a_{mn} \cdot x_n & = b_m & (m) \end{array}\right.$

Прямой ход:

$\begin{array}{ccc} (2)\to (2)-(1) \cdot ( \frac {a_{21}}{a_{11}}) &:& a_{22}^{\prime} \cdot x_2 + a_{23}^{\prime} \cdot x_3 + \ldots + a_{2n}^{\prime} \cdot x_n = b_2^{\prime} \\ (3)\to (3)-(1) \cdot ( \frac {a_{31}}{a_{11}}) &:& a_{32}^{\prime} \cdot x_2 + a_{33}^{\prime} \cdot x_3 + \ldots + a_{3n}^{\prime} \cdot x_n = b_3^{\prime} \\ \ldots & & \\ (m)\to (m)-(1) \cdot ( \frac {a_{m1}}{a_{11}}) &:& a_{m2}^{\prime} \cdot x_2 + a_{m3}^{\prime} \cdot x_3 + \ldots + a_{mn}^{\prime} \cdot x_n = b_n^{\prime} \\ (3)\to (3)-(2) \cdot ( \frac {a_{32}^{\prime}}{a_{22}^{\prime}}) &:& a_{33}^{\prime\prime} \cdot x_3 + \ldots + a_{3n}^{\prime\prime} \cdot x_n = b_3^{\prime\prime} \\ \ldots & & \\ (m)\to (m)-(m-1) \cdot ( \frac {a_{m,n-1}^{(m-2)}}{a_{m-1,n-1}^{(m-2)}}) &:& a_{mm}^{(m-1)} \cdot x_m + \ldots + a_{mn}^{(m-1)} \cdot x_n = b_m^{(m-1)} \end{array}$

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

Пример

Покажем, как методом Гаусса можно решить следующую систему:

$\left\{\begin{array}{ccc}2x + y - z &=& 8 \\ -3x - y + 2z &=& -11 \\ -2x + y + 2z &=& -3 \end{array}\right.$

Обнулим коэффициенты при $x\!$ во второй и третьей строчках. Для этого вычтем из них первую строчку, умноженную на $\textstyle-\frac{3}{2}\!$ и $-1\!$ , соответственно:

$\left\{\begin{array}{rcc} 2x + y - z &=& 8 \\ \frac{1}{2}y + \frac{1}{2}z &=& 1 \\ 2y + z &=& 5 \end{array}\right.$

Теперь обнулим коэффициент при $y\!$ в третьей строке, вычтя из неё вторую строку, умноженную на $4\!$ :

$\left\{\begin{array}{rcc} 2x + y - z &=& 8 \\ \frac{1}{2} y + \frac{1}{2} z &=& 1 \\ -z &=& 1 \\ \end{array}\right.$

В результате мы привели исходную систему к треугольному виду, тем самым закончив первый этап алгоритма.

На втором этапе разрешим полученные уравнения в обратном порядке. Имеем:

$z = -1\!$ из третьего;

$y = 3\!$ из второго, подставив полученное $z\!$

$x = 2\!$ из первого, подставив полученные $z\!$ и $y\!$ .

Таким образом исходная система решена.

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