Моделирование электротехнических систем и систем автоматики

где xk–1 − значение корня на предыдущей итерации; xk − значение корня на текущей итерации.

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

Рис. 16. Графическое изображение метода поиска корней системы двух линейных уравнений

На рис. 16 изображены две линейные алгебраические зависимости, решение которых нужно найти приближенным итерационным методом. Для этого произвольно выберем начальное значение корня уравнения (1). Пусть это будет аргумент Х0. Вставим это значение в уравнение (2). В этом случае первое приближенное решение этой системы будет Х0, Х1 при значении обеих функций, равном Y0. При этом решении системы интервал от Х0 до Х1 будет намного превышать заданный диапазон точности решения, поэтому на следующем этапе поиска решения вставим значение аргумента Х1 в уравнение (2) и проецируя значение Y1 на линию уравнения (1) получим при данном значении функции решение системы на следующей итерации равное Х1, Х2. Диапазон между корнями Х1 и Х2 при этом значительно сужается.

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

41

Реализация этого принципа определения на ЭВМ корней системы алгебраических уравнений представлена алгоритмом

(рис. 17).

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

42

После запуска этого алгоритма вводятся следующие параметры:

n − число уравнений в системе (порядок системы);

A(n, n + 1) − двумерный массив коэффициентов системы; X(n) − массив начальных значений корней системы;

ε − заданная точность решения на участке приращения аргумента ∆x.

После ввода исходных данных открывается вложенный цикл, с заданным числом повторений по параметрам K и J, которые изменяются в пределе от 1 до n.

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

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

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

6.3. Принцип и алгоритм численного дифференцирования

Дифференциалом dy функции f(x) называется приращение ординаты касательной этой функции, проведенной из точки начала этого приращения (рис. 18).

43

Рис. 18. Изменение дифференциала dy и приращения функции ∆y на участке приращения аргумента ∆х

Производной функции называют отношение ее дифференциала dy (приращения ординаты касательной) к приращению аргумента dx:

Геометрически производная – это тангенс угла α наклона касательной к приращению аргумента ∆x (см. рис. 18).

Численное дифференцирование функции производится через вычисление полинома, структура которого зависит от числа точек (узловых значений функции) на участке приращения аргумента. Эти значения функции на участке дифференцирования называют «узлами».

Если число этих узлов равно трем, то вычисление полинома производится по схеме, представленной на рис. 19, а.

Рис. 19. Схема вычисления полинома при численном дифференцировании функции

44

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

где x0 − начальное значение аргумента при дифференцировании; x − конечное значение аргумента при дифференцировании; h − шаг дифференцирования.

Сучетом выражения (6) при дифференцировании функции

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

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

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

Алгоритм численного дифференцирования функции показан на рис. 20. После пуска этого алгоритма необходимо задать следующие параметры: число точек (узлов) дифференцирования N, шаг дифференцирования функции h и диапазон ее диф-

Рис. 20. Алгоритм численного дифференцирования функции

46

ференцирования от x0 до x. На основе этих данных по формуле (6) вычисляется параметр p, который затем вставляется в формулу конкретного полинома при вычислении дифференциала. Далее алгоритм анализирует величину параметра N и в зависимости от его значения требует ввода соответствующего числа значений функций в выбранных точках. После ввода необходимого числа значений функций алгоритм по соответствующей формуле полинома вычисляет величину дифференциала.

6.4. Принцип и алгоритм численного вычисления определенного интеграла

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

трапеций.

По методу «прямоугольника» (рис. 21, а) площадь подинтегральной фигуры заменяется суммой площадей прямоугольников, основание которых равно шагу интегрирования Н, а высота определяется значением функции в точке очередного приращения аргумента.

Величина шага приращения аргумента определяется как:

где В – верхний предел интегрирования; А – нижний предел интегрирования; N – число участков от А до В.

47

Рис. 21. Схема вычисления определенного интеграла (а) методом «прямоугольника»; (б) методом трапеций

Элементарная площадь прямоугольника подинтегральной функции

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

N

S = ∑[ f (A) + f (A + H ) +... + f (A + H (N −1)]H. (28)

i=1

Алгоритм вычисления определенного интеграла методом прямоугольника представлен на рис. 22.

После пуска этого алгоритма необходимо задать пределы интегрирования функции (параметры А и В) и число N интервалов разбиения диапазона интегрирования функции. Параметром ε вводится заданная точность вычисления определенного интеграла.

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

48

Рис. 22. Алгоритм интегрирования функции методом «прямоугольника»

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

49

суммарной площади элементарных прямоугольников подинтегральной функции в пределах интервала ее интегрирования.

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

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

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

i=1

Алгоритм для вычисления определенного интеграла методом «трапеций» представлен на рис. 23.

Вотличие от предыдущего метода этот алгоритм сразу

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

50