П.Г. Колпахчьян, И. Б. Подберезная, С.В.Чамлай, Д.В. Батищев

П.Г. Колпахчьян, И. Б. Подберезная, С.В.Чамлай, Д.В. Батищев

 

Методы вычислений в задачах электроаппаратостроения, электрооборудования и электрического транспорта: Методические указания к проведению вычислительной практики. Юж.-Рос. гос. техн. ун-т(НПИ). Новочеркасск: ЮРГТУ(НПИ), 2012. - 75с.

 

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

Указания предназначены для студентов первого и второго курсов дневной и заочной форм обучения направления 14040062 «Электроэнергетика и электротехника» ФГС ВПО (III поколения) профилей №9 Электрические и электронные аппараты, №12 Электрический транспорт, №16 Электрооборудование и электрохозяйство предприятий, организаций и учреждений.

Методические указания утверждены на заседании кафедры «Электрические и электронные аппараты», «Электрический транспорт» протокол №7, от «29»февраля 2012г.

 

 

Ó Южно-Российский государственный

технический университет, 2012

Ó Колпахчьян П.Г., Подберезная И. Б.,

Чамлай С.В., Батищев Д.В., 2012

Оглавление

Предисловие. 5

Глава 1. Решение алгебраических уравнений и их систем.. 7

1.1. Решение систем линейных алгебраических уравнений. 7

1.1.1. Итерационные методы решения систем линейных алгебраических уравнений 7

1.1.2. Факторизация и типовые схемы решений. 9

1.1.3. Метод Гаусса и LU — разложение. 11

1.1.4. Метод Гаусса-Жордана обращения матриц. 13

1.1.5. Метод квадратного корня (Холесского) 13

1.1.6. Метод вращений. 15

1.1.7. Итерационное уточнение. 16

1.1.8. Решение преопределенных систем линейных алгебраических уравнений 17

1.2 Решение систем нелинейных алгебраических уравнений. 18

1.2.1.Метод последовательных приближений. 19

1.2.2. Метод Ньютона. 20

1.2.3. Метод Ньютона по параметру. 21

Глава 2. Интерполяция зависимостей. 23

2.1. Интерполяция каноническим полиномом. 23

2.2. Интерполяция полиномом Лагранжа. 25

2.3. Интерполяция полиномом Ньютона. 26

2.4. Применение интерполяции для решения уравнений. 28

2.5. Интерполяция сплайнами. 30

Глава 3. Метод наименьших квадратов. 34

3.1. Общие положения. 34

3.2. Степенной базис. 36

3.3. Базис в виде классических ортогональных полиномов. 36

3.4. Базис в виде ортогональных полиномов дискретной переменной. 38

3.5. Линейный вариант метода наименьших квадратов. 39

3.6. Сглаживание экспериментальных данных с ошибками. 40

Глава 4. Определенные интегралы.. 42

4.1. Классификация методов. 42

4.2. Методы прямоугольников. 43

4.3. Апостериорные оценки погрешностей по Рунге и Эйткену. 46

4.4. Метод трапеций. 48

4.5. Метод Симпсона. 49

4.6. Методы Ньютона-Котеса. 51

4.7. Вычисление интегралов с заданной точностью.. 52

4.8. Применение сплайнов для численного интегрирования. 54

4.9. Методы наивысшей алгебраической точности. 54

4.10 Несобственные интегралы.. 58

4.11 Вычисление кратных интегралов. 59

Глава 5. Решение обыкновенных дифференциальных уравнений. 60

5.1. Типы задач для обыкновенных дифференциальных уравнений. 60

5.2. Метод Эйлера. 61

5.3. Методы Рунге-Кутта. 63

5.4. Метод Рунгe-Кутта-Мерсона. 66

5.5. Методы Адамса-Башфорта и Адамса-Маултона. 67

5.6. Методы Гира. 70

Библиографический список…………………………………………………….73

 

 

 

 

Предисловие

Настоящее пособие предназначено для студентов первого, второ­го курсов направления 14040062 Электроэнергетика и электротехника при прохождении вычислительной практики. От­дельные разделы могут быть использованы в курсах «Математическое моделирование», «Переходные режимы работы электрооборудования» и т.д.

Продолжительность вычислительной практики две недели. Прово­дится она в вычислительных центрах и вычислительных лабораториях. Распределение студентов на практику оформляется приказом по университету не позднее, чем за месяц до ее начала. Перед началом практи­ки студенты должны получить на кафедре бланки командировочных удостоверений, пройти необходимый инструктаж.

Цели и задачи вычислительной практики Основная цель вычислительной практики — это закрепление теоретических знаний, по­лученных в университете, и овладение практическими навыками реше­ния конкретных задач на ЭВМ с использованием численных методов. Теоретическую основу практики составляют курсы "Информатика", "Алгоритмические языки и программирование", "Высшая математика". В задачи практики вхо­дит:

- изучение методов решения систем линейных и нелинейных урав­нений, систем дифференциальных уравнений первого порядка, а также знакомство с методом наименьших квадратов;

- изучение алгоритмического языка PASCAL;

- составление и отладка программ для перечисленных методов на алгоритмическом языке PASCAL.

Содержание практики. В период вычислительной практики каждый студент должен:

- пройти на кафедре и на месте проведения практики инструктаж
по технике безопасности;

- ознакомиться с литературой по методам решения указанных вы­ше задач (см. библиографический список);

- изучить алгоритмический язык PASCAL для работы на компьютерах;

- составить и отладить на алгоритмическом языке PASCAL программы, реализующие:

1) методы решения систем обыкновенных дифференциаль­ных уравнений первого порядка с заданными начальными значениями;

2) методы решения систем линейных уравнений общего вида;

3) методы решения систем нелинейных уравнений;

4) аппроксимацию функции по методу наименьших квадратов, интерполяцию и сглаживание;

5) вычисление функций с помощью разложения в степенные ряды;

6) расчет переходных процессов в электрических цепях;

- оценить точность полученных результатов, составить отчет по практике.

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

Отчет должен со­держать следующие разделы:

- математическая сущность используемых методов;

- алгоритмы методов, указанных в индивидуальном задании сту­дента; пояснения к блок-схемам программ и описание программ, реализующих эти методы;

- индивидуальный вариант расчета;

- выводы (сравнение результатов, полученных различными мето­дами; затраты машинного времени; простота алгоритма; объем используемой памяти ЭВМ);

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

Подведение итогов практики. В период вычислительной практики для всех студентов является обязательной ежедневная работа в вычислительной лаборатории. Студенты должны вести дневник о ежедневно выполненной работе, по которому составляется отчет. Составленный студентом отчет вместе с дневником и другими документами должен быть представлен на кафедру для проверки руководителем. После про­верки материалов практики студент защищает отчет.

Глава 1

Метод вращений

Матрицы вращений позволяют реализовать упорядоченное исклю­чение переменных.

Построим ортогональную матрицу вращений R₂₁ так, чтобы при левом умножении матрицы, обратной к ней, на матрицу А она обра­щала в ноль элемент агь (исключение переменной x1 из второго урав­нения) затем матрицу R₃₁ для обращения в ноль а и так далее, пока не будут обращены в ноль поддиагональные элементы первого столбца матрицы А:

……………………………………….

Выполним подобную последовательность операций для всех столбцов матрицы А. Получим факторизацию, на которой основан ме­тод вращений (Гивенса):

 

Схема выполнения операций выглядит следующим образом

 

Очевидно, что в результате обратится в ноль элемент .

Если необходимо решить одну систему уравнений с одной правой частью, то можно использовать следующий алгоритм решения:

1) Для k=1,2,3,…, n -1 выполнить пп. 2,3,4,5.

2) Для i=k+1, k+2,…,n выполнить пп. 3,4,5.

3) Вычислить элементы c и s матрицы вращения

4) Умножить матрицу A слева на , для чего

a) положить

b) для j=k+1,k+2,…,n выполнить

5) Умножить Матрицу на вектор правых частей, для чего выполнить

6) Решить полученную систему уравнений в верхней треугольной форме методом обратной подстановки.

7) Конце алгоритма. Матрица A содержит в верхнем треугольнике матрицу U, вектор B – пересчитан.

 

Каждый шаг исключений переменной методом вращений из очередного k-го столбца системы уравнений состоит из нескольких малых шагов, на каждом из которых необходимо строит матрицу .

Следует обратить внимание на то, что в отличии от метода Гаусса пересчету на кадом малом шаге подлежат так же и элементы ведущей строки c индексами столбцов j, i, k.

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

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

 

Итерационное уточнение

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

1) Решить исходную систему уравнений .

2) Вычислить вектор невязок . Если все , где - желаемая точность, то расчет закончен и есть решение, иначе выполнить пп. 3,4,5.

3) Решить систему уравнений относительно - вектора итерационного уточнения.

4) Вычислить уточненное решение .

5) Повотрить пп. 2,3,4,5.

6) Конец алгоритма.

 

Приведенный алгоритм вытекает из следующих соотношений

.

Т.к. AX=B, имеем:AX=R.

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

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

 

Метод Ньютона

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

Рассмотрим простой пример (рис. 1.1).

Рис. 1.1. Итерация метода Ньютона для

 

Поскольку где — начальное приближение, то и

dw(x₀)dx можно получить новое приближение .

Продолжая итерационный процесс, можно с требуемой точно­стью риблизиться к одному из решений, например, = 1,73205...

Расчетная формула для метода Ньютона может быть получена, ес­ли представить w(x) в окрестности текущего приближения х_k в виде ряда Тейлора

 

 

и ограничиться линейными членами, тогда в матричной форме получим

 

 

где

Применительно к СНАУ W(X) = 0 получим следующий алгоритм:

1) Выбрать начальный вектор положить к = 0, =

2) Вычислить вектор . Если все < , где — за­данная точность расчета, то получено решение, расчет окончен..

Если к > 1 и max , то итерационный

процесс расходится, расчет завершить аварийно.

3) Построить матрицу Якоби

 

и вычислить значения всех производных в точке

4) Решить систему уравнений, определив вектор поправок АХ^:

5) Вычислить новое приближение

и положить к = к + 1.

6) Если к > к_тах, где k_max — заданное предельное число итераций, то аварийно завершить расчет, иначе перейти к п. 2.

7) Конец алгоритма.

Метод Ньютона при начальном приближении, близком к некото­рому решению, часто обладает устойчивой квадратичной сходимостью. При плохой начальной точке имеет место расходящийся итераци­онный процесс. Метод Ньютона расходится, если матрица Якоби пло­хо обусловлена в окрестности решения. Часто перед использованием метода Ньютона выполняют несколько итераций, например методом последовательных приближений, для того чтобы иметь "хорошее" на­чальное приближение.

В качестве косвенного критерия расходимости итерационного про­цесса можно использовать изменение знака якобиана — определите­ля матрицы Якоби. Однако это условие будучи достаточным, не явля­ется необходимым. Якобиан может быть вычислен, как побочный про­дукт решения методом Гаусса системы из п. 3 алгоритма.

Метод Ньютона по параметру

Метод Ньютона по параметру относится к классу квазиньюто­новских методов и предусматривает расчет нового исходного прибли­жения по формуле

где — параметр, выбираемый на каждой итерации.

При — 1 метод совпадает с обычным методом Ньютона.

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

Метод Бройдена предусматривает выбор параметра в про­цессе поиска по направлению вектора поправок минимума од­ной из норм вектора невязок, например, евклидовой:

В простейшем случае можно попытаться выполнить последова­тельность расчетов Ф(Х), уменьшая шаг вдвое, т.е. полагая 0.5,0-25,... пока Ф(Х) не уменьшится, или не станет малой ве­личиной. Более сложные варианты предусматривают аппроксимацию сечения по направлению поверхности функции невязок выпук­лой квадратичной функцией и выбор , обеспечивающего минимум этой функции.

Вычислим значение Ф₀ и, выполнив предварительно шаг метода Ньютона, вычислим Если , то поиск на отрезке [0,1] начи­нается с а = 1, в противном случае с а = 0. Зададимся некоторым и вычислим значения функции Ф в точках (1 — ) в первом случае и в точках во втором случае, где k = 1, 2,... Поиск прекращается при выполнении условия

В окрестности k-й точки производится квадратичная аппроксима­ция функции Ф(Х) и определяется из условия ее минимума (равенство нулю производной) значение , равное

 

где — точки, в которых известно значение аппроксими­рующей функции

Метод Матвеева предусматривает выбор оптимального парамет­ра следующим образом:

 

 

Значение определяется, например, по следующей формуле:

 

 

где — i -e решаемое уравнение на к -й итерации, а второй сомножитель представляет собой максимальный по модулю элемент вектора-столбца, полученного умножением матрицы вторых производ­ных системы W(X) на вектор поправок.

Для метода Матвеева необходимо, чтобы функции Wi были дважды дифференцируемыми по

Метод Матвеева требует большего по сравнению с методом Нью­тона объема вычислений и памяти, однако обладает устойчивой сходи­мостью при не очень хороших исходных приближениях. Доказано, что данный метод сколь угодно приближается к решению, даже если яко­биан стремится к нулю, т.е. матрица Якоби очень плохо обусловлена в окрестности решения.

Глава 2

Интерполяция зависимостей

Интерполяция сплайнами

Полиномиальная интерполяция не всегда дает удовлетворитель­ные результаты при аппроксимации зависимостей. Так, например, при представлении полиномами резонансных кривых колебательных си­стем большая погрешность возникает на концах ("крыльях") этих кри­вых. Несмотря на выполнение условий Лагранжа в узлах, интерполя­ционная функция может иметь значительное отклонение от аппрокси­мируемой кривой между узлами. При этом повышение степени интер­поляционного полинома приводит не к уменьшению, а к увеличению погрешности. Возникает так называемое явление волнистости [4].

Для проведения гладких кривых через узловые значения функции чертежники используют упругую металлическую линейку, совмещая ее с узловыми точками. Математическая теория подобной аппроксима­ции развита за последние тридцать лет [5, 6] и называется теорией сплайн-функций. Разработано и обширное программное обеспечение для применения сплайнов в различных областях науки и техники [7, 1].

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

(2.18)

Функцию и будем использовать для аппроксимации зависимости f(x), заданной в узлах ,...,х_п значениями .

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

(2.19)

где и — коэффициенты сплайна, определяемые из дополни­тельных условий; i = 1,2,..., п — номер сплайна.

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

Коэффициенты сплайнов определяются из условий сшивания со­седних сплайнов в узловых точках:

1) Равенство значений сплайнов и аппроксимируемой функ­ции f(x) в узлах — условия Лагранжа

(2.20)

2) Непрерывность первой и второй производных от сплайнов в уз­лах

(2.21)

(2.22)

Кроме перечисленных условий необходимо задать условия на кон­цах, т.е. в точках и. В общем случае эти условия зависят от кон­кретной задачи. Часто используют условия свободных концов сплай­нов. Если линейка не закреплена в точках вне интервала [ , ], то там она описывается уравнением прямой, т.е. полиномом первой степени. Следовательно, исходя из условий (2.22), т.е. непрерывности вторых производных сплайнов на концах интервала, запишем соотношения:

(2.23)

- (2.24)

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

Получим алгоритм определения коэффициентов кубических сплай­нов из условий (2.20 — 2.24). Условия (2.20) в узлах x_i-1 и x _i после подстановки i -гo сплайна принимают вид

a_i = ƒ_{i-1 ,} (2.25)

a_i + b_i h_i + c_ih_i² + d_ih_i³ = ƒ_{i ,}(2.26)

где h_i = x_i — x_{i- 1}, 1 ≤ i ≤ n.

Продифференцируем дважды сплайн (2.19) по переменной x:

φ′_i(x) = b_i + 2c_i (x— x_i-1) + 3d_i (x— x_i-1)², (2.27)

φ″_i(x) = 2c_i + 6d_i (x— x_i-1). (2.28)

Из условий непрерывности производных (2.21) и (2.22) при пере­ходе в точке x_iот i-го к (i₁)сплайну с учетом выражений (2.27) и (2.28) получим следующие соотношения:

b_i + 2c_i h_i + 3d_ih_i² = b_{i+1 ,}(2.29)

c_i + 3d_ih_i = c_i + 1. (2.30)

И, наконец, из граничных условий (2.23) и (2.24) на основании вы­ражения для второй производной (2.28) получим, что

c₁ = 0, (2.31)

c_n + 3d_n h_n = 0. (2.32)

Соотношения (2.25), (2.26), (2.29) — (2.32) представляют собой полную систему линейных алгебраических уравнений относительно ко­эффициентов сплайнов a_i, b_i, c_i и d_i. Но прежде чем решать эту систе­му, выгодно преобразовать ее так, чтобы неизвестными была только одна группа коэффициентов c_i.

Из уравнений (2.30) коэффициенты d_iвыразим через коэффициен­ты c_i:

c_i+1 —c _i

d_i = ————. (2.33)

3 h_i

Объединяя уравнения (2.25) и (2.26) с уравнением (2.33), предста­вим коэффициенты b_iтакже через коэффициенты c_i:

f _i — f _i-1(c _i+1 + 2c _i)h

b_i = ———— - —————. (2.34)

h_i 3

После подстановки выражений (2.33) и (2.34) в соотношение (2.29) получим уравнение, в которое входят только неизвестные коэффици­енты c_i. Для симметричности записи в полученном уравнении умень­шим значение индекса i на единицу:

h_i-1 c_i-1+2(h_i-1+h_i)c_i + h_ic_i+1 = 3 (2.35)

где 2≤ i ≤ n.

При i = n, учитывая условие свободного конца сплайна, в уравне­ний (2.35) следует положить

c_n+1 = 0. (2.36)

Таким образом, п — 1 уравнение вида (2.35) вместе с условиями (2.31) и (2.36) образует систему линейных алгебраических уравнений для определения коэффициентов c_i. Коэффициенты d _i и b _i вычисляют­ся после нахождения c_i по формулам (2.33) и (2.34), коэффициенты a_i равны значениям аппроксимируемой функции в узлах в соответствии с формулой (2.25).

В каждое из уравнений типа (2.35) входит только три неизвестных с последовательными значениями индексов c_i-1, c_i, c_i+1. Следователь­но, матрица системы линейных алгебраических уравнений относитель­но c_i является трехдиагональной, т.е. имеет отличные от нуля элементы только на главной и двух примыкающих к ней диагоналях. Для реше­ния систем с трехдиагональной матрицей наиболее эффективно при­менять так называемый метод прогонки, являющийся частным случаем метода Гаусса.

Рассмотрим алгоритм метода прогонки применительно к решаемой задаче. Для сокращения записи уравнение (2.35) представим в виде

h_i-1 c_i-1+s_i c_i + h_ic_i+1 = r_i, (2.37)

где s_i = 2 (h_i-1 - h_i ),

r_i = 3 . (2.38)

Так же, как и метод Гаусса, метод прогонки разделяется на два эта­па — прямой и обратный ходы. В процессе прямого хода метода про­гонки вычисляют значения коэффициентов линейной связи каждого предыдущего неизвестного c_i с последующим c_i+1.

Из уравнения (2.37) при i = 2 с учетом граничного условия (2.31) установим связь коэффициента с₂ с коэффициентом с₃:

 

c₂ = k₂ - l₂с_3, (2.39)

где k₂ и l₂ — прогоночные коэффициенты,

 

. r ₂ h₂

k₂ = —, l₂ = —.

S ₂ S₂

Затем, подставляя выражение (2.39) в уравнение (2.37) при i = 3, получим линейную связь коэффициента с₃ с коэффициентом с₄:

c₃ = k₃ - l₃с_4,

Поступая аналогичным образом для любых соседних коэффициентов с номерами i и i + 1, можно установить линейную связь в виде

c_i = k_i - l_iс_i+1 (2.40)

В процессе выполнения прямого хода метода прогонки необходимо вы­числить значения всех прогоночных коэффициентов k _i и l_i, для кото­рых получим рекуррентные соотношения. Подставим формулу связи

(i -1)-го и i -го коэффициентов

c_i-1 = k_i-1 - l_i-1с_i

в уравнение (2.37), в результате получим

r_i - h_i-1k_i-1h_i

c_i = ――――― - ————— c_i+1.

s_i - h_i-1l_i-1 s_i - h_i-1 l_i-1

Сравнение последнего соотношения с формулой (2.40) позволяет за­писать рекуррентные соотношения для определения прогоночных ко­эффициентов k_i и l_i:

r_i - h_i-1k_i-1h_i

k_i = ――――――, l_i = ——————. (2.41)

s_i - h_i-1l_i-1 s_i - h_i-1 - l_i-1

 

Учитывая граничное условие (2.31) и соотношение

c₁ = k₁ - l₁с_2,

а также полагая c₂ ≠0, задаем начальные коэффициенты k₁ = 0 и l₁ = 0. Затем по формуле (2.41) вычислим все п пар прогоночных ко­эффициентов k _i и l_i..

На основании соотношения c_n = k_n - l _nc_n+1 и граничного условия (2.36) получим, что

c_n = k_n. (2.42)

далее последовательно применим формулу (2.40) при i = n — 1, n — 2,...., 2 и вычислим значения искомых величин c_n-1, c_n-2,…, c_{2 .}

Эта процедура называется обратным ходом метода прогонки.

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

После определения всех коэффициентов с_i другие коэффициенты сплайна вычисляются по формулам (2.25), (2.33) и (2.34). Аппрокси­мирующую функцию φ(x) можно рассчитать с помощью соотношения (2.19) в любой точке x на интервале [ х₀, х_n ].

Глава 3

Метод наименьших квадратов

Общие положения

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

Обозначим узлы исходной таблицы данных через x _j, где 0 ≤ j ≤ n — номер узла. Считаем известными значения эксперимен­тальных данных в узловых точках ƒ(x_j) = ƒ_j. Введем непрерывную функцию φ (х) для аппроксимации дискретной зависимости f(x _j). В узлах функции φ (х) и f (x)будут отличаться на величину ε _j = φ (x _j) — - f(x_j). Отклонения ε _j могут принимать положительные и отрица­тельные значения. Чтобы не учитывать знаки, возведем каждое откло­нение в квадрат и просуммируем квадраты отклонений по всем узлам:

n n

Q = ∑ ε _j ² = ∑[φ(x _j) - ƒ(x _j)]² (3.1)

j=0 j=0

Метод построения аппроксимирующей функции φ (х)из условия минимума величины Qполучил название метода наименьших квадра­тов (МНК).

Наиболее распространен способ выбора функции φ (х)в виде ли­нейной комбинации

φ (х) = c₀ φ ₀ (x)+ c₁ φ ₁ (x) +... + c_m φ _m(x), (3.2)

φ ₀ (х), φ ₁ (x),..., φ _m (х) — базисные функции; m ≤ n;

c₀(x), c₁ (x),.. c_m (x)- коэффициенты, определяемые при минимизации величины Q.

Математически условия минимума квадратов отклонений Q запи­шем, приравнивая к нулю частные производные от Qпо коэффициентам 0 ≤ k ≤ m:

n

∂ Q⁄∂ c₀ = 2∑ [c₀ φ ₀ (x _j) + c₁ φ ₁(x _j) +…+

j=0

+ c_m φ _m (x _j) - ƒ_j] φ ₀ (x _j) = 0,

n

∂ Q⁄∂ c₁ = 2∑ [c₀ φ ₀ (x_j)+c₁φ₁(x_j) +…+ (3.3 )

j=0

+ c_m φ _m (x _j) - ƒ_j] φ ₁ (x _j) = 0,

…………………………………….

n

∂ Q⁄∂ c_m = 2∑ [c₀ φ ₀ (x _j) + c₁ φ ₁(x _j) +…+

j=0

+ c_m φ _m (x _j) - ƒ_j] φ _m (x _j) = 0,

Из системы линейных алгебраических уравнений (3.3) определя­ются все коэффициенты с‌‍_k. Система (3.3) называется системой нор­мальных уравнений. Матрица этой системы имеет следующий вид:

(3.4)

и называется матрицей Грамма. Элементы матрицы Грамма являются скалярными произведениями базисных функций:

n

(φ _j, φ _k) = ∑ φ _j (x _j) φ _k(x _j). (3.5)

j=0

 

Расширенная матрица системы уравнений (3.5) получится добав­лением справа к матрице Грамма столбца свободных членов

(3.6)

,

где скалярные произведения, являющиеся элементами столбца, опре­деляются аналогично (3.5):

(3.7)

n

(φ _j, ƒ) = ∑ φ _j (x _j) ƒ _j

j=0

Отметим основные особенности матрицы Грамма, полезные при про­граммной реализации алгоритмов МНК:

1) Матрица симметрична, т.е. a‍‍ _{i‍‍‍‍‍‍ j} = a‍‍_{‍‍‍ j i}, что позволяет сократить объем вычислений при заполнении матрицы.

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

3) Определитель матрицы будет отличен от нуля, если в качестве базиса выбраны линейно независимые функции φ‌_k (x), при этом система (3.3) имеет единственное решение

При обработке экспериментальных данных, определенных с погреш­ностью е в каждой узловой точке, обычно начинают с аппроксима­ции функцией φ‌_k (x), представимой одной-двумя базисными функци­ями. После определения коэффициентов c‍_kвычисляют величину Q по формуле(3.1). Если получится, что √Q > e, то необходимо расширить базис добавлением новых функций φ‌_k (x).Расширение базиса необхо­димо осуществлять до тех пор, пока не выполнится условие √Q ≈ e.

Выбор конкретных базисных функций зависит от свойств аппрок­симируемой функции ƒ(x), таких как периодичность, экспоненциаль­ный или логарифмический характер, свойства