Файл: Численные методы и математическое моделирование.pdf

13 данной функции. Следовательно, надо найти экстремум и потом искать корни в двух интервалах – от -∞ до точки экстремума и от точки экстре- мума до ∞. Данный подход может быть распространён и на уравнения с большим количеством корней.
Однако часть корней могут быть иррациональными числами (ком- плексными). Чтобы выполнить данную проверку, надо сравнить значе- ния функции на концах отрез- ков, которые содержат корни: если эти числа имеют разные знаки, то корень существует.
Как это выглядит графически, показано на примере кубиче- ского уравнения (которое мо- жет иметь до 3 корней) (рис. 3).
Как видно, области первого и второго корней лежат выше нуля оси Х. Таким образом, они отсутствуют среди рациональ- ных чисел.
Как говорилось во введении существуют два подхода в реализации этих решений: проекционный и конечно-разностный. Рассмотрим их оба на примере нахождения корня уравнения.
В проекционных методах строятся упрощённые зависимости, кото- рые заменяют исследуемую функцию и позволяет найти новое прибли- жение к решению задачи. Данный подход реализован в методах итера- ций и касательных. Рассмотрим эти метод более подробно.
Алгоритмы и их реализация в Excel
Метод касательных (хорд)
Метод является проекционным и в качестве проекции используется производная от функции. Функция в заданной точке рассматривается как линейная, соответствующая производной, которая является каса- тельной к решаемому уравнению в этой точке и задаётся тангенсом угла наклона касательной. Зная координаты точки и угол наклона касатель-
Y
X
0
X
к3
X
э2
X
э1
X
к1
X
к2
X
к3
Рис. 3. Существующий и несуществующий корни

14 ной строим линейное уравнение и для него находим корень (точка пе- ресечения касательной с осью Х), который является следующим приближением к решению.
Графическая схема метода показана на рис. 4. Касательные можно заменять хордами, прове- дёнными через две точки, равно- отстоящие от заданной точки.
Метод очень чувствителен к вы- бору начальной точки.
Алгоритм метода показан на рис.
5. Рассмотрим его реализацию:
Обычно имеются два вида уравне- ний для решения:
 уравнение с известным алгебраи- ческим выражением для произ- водной, при этом задача сводится лишь к нахождению координат точки пересечения производной с осью Х (об этом далее);
 алгебраического выражения про- изводной нет, то строим прямую одним из двух методов: можно использовать программу числен- ного вычисления производной в заданной точке или воспользо- ваться методом хорд, когда через две точки строится прямая и в точке её пересечения с осью Х
Y
X
0
X
1
X
0
X
2
Рис.4 Графическая реализация ме- тода касательных (хорд)
Рис.5. Алгоритм метода касательных (хорд)

15 возникает следующая точка приближения к корню. Далее проце- дура повторяется. В этом случае надо на первом шаге задать не одну, а две точки или организовать поиск второй точки по какой- либо процедуре, например, с заданным шагом вперёд по оси Х.
Чтобы реализовать эти методы вспомним процедуру нахождения не- известных коэффициентов линейного уравнения вида Y=a+b·X на ос- новании имеющихся данных:
 при известном угле наклона b (значение нашей производной) и на основании известных координат заданной точки (X
0
, Y
0
) находим a по формуле:
b
X
Y
a



0 0
;
 по двум точкам (X
1
, Y
1
) и (X
2
, Y
2
) находим
1 2
1 2
X
X
Y
Y
b



и потом
b
X
Y
a



, используя одну из имеющихся точек.
Из полученного уравнения прямой находим точку его пересечения с осью Х по условию Y=0, что приводит к формуле:
b
a
X


. Это зна- чение является приращением (Х) к предыдущему значения Х
i-1
для вы- числения следующего приближением к корню Х
i
= Х
i-1
+ Х. Повторяем процедуру, до тех пор, пока разница между двумя последними Х
i
не ста- нет меньше заданной погрешности erf.
Рассмотрим реализацию метода в пакете электронных таблиц MS
Excel. Открываем таблицу и начинаем готовить решение. В качестве уравнения используем полином второй степени, для этого сначала на листе создаём таб- личку с коэффици- ентами полинома и необходимы ком- ментарии к реше- нию (рис. 6).
Построим таб- лицу функции для подбора коэффициентов полинома по условию, что один из корней лежит в удобной для наблюде- ния на графике области. Сначала готовим за- головок и шапку таблицы (рис. 7), потом в столбце «Х» записываем два числа (начальное значение и следующее для задания шага табу- ляции), выделяем их мышкой и растягиваем за
Рис. 6. Таблица коэффициентов
Рис. 7 Заголовок и шапка таблицы

16 чёрный маркёр в правом нижнем углу выделе- ния черным крестиком курсора мышки до жела- емого конечного значения (рис. 8).
В столбец «Y=F(X)» записываем формулу нашей функции:
=$B$2+$B$3*A8+$B$4*A8^2
, где
$B$2
,
$B$3
,
$B$4
– коэффициенты полинома, ко- торые берутся из строго определённых ячеек
(см. рис. 6), поэтому закрепляем их абсолют- ными ссылками (знаки
$
ставится с помощью одного щелчка клавиши
F4
);
A8
– ссылка на входной параметр «X» явля- ется относительной ссылкой. Растягиваем формулу до конца таблицы используя чёрный хвостик в правом нижнем углу рамки вокруг ячейки с формулой, который хватаем черным крестиком курсора мышки, и по- лучаем таблицу для построения графика (рис. 9а). По её данным строим график (рис. 9б). а. б.
Рис. 9. График функции
При необходимости можно, изменяя коэффициенты полинома, до- биться, чтобы один из корней уравнения лежал в нашем диапазоне по
«Х», например, в районе 3,5 как видно из графика.
В нашем случаем известны и уравне- ние и его производная, поэтому решение реализуется просто. Готовим заготовку таблицы решения (рис. 10) и заполняем
Рис. 8. Построение прогрессии
Рис. 10. Заголовок и шапка расчётной таблицы

17 первую строку данными. В столбец «Х» ставим начальное приближение для поиска корня и заливаем эту ячейку цветом для напоминания, что здесь надо ввести данные. Столбец «Y=F(X)» содержит ту же формулу, что и первая таблица для табулирования функции. Её можно просто скопировать из первой таблицы (
Ctrl+C
) и вставить во вторую (
Ctrl+V
).
Вычисляем производную для полинома, формула которой имеет сле- дующую запись в таблице Excel:
=$B$3+2*$B$4*D8
. Остаётся внести фор- мулу для определения следующего шага по «Х» во второй строке таблицы. Фор- мула получается из треугольника, кото- рый строится от заданной начальной точки и вдоль касательной до пересече- ния с осью «Х» (рис. 11). Отрезок (Х
0
, Х
1
) прибавляется к Х
0
с учётом знака. Дан- ная точка (Х
1
) будет следующим прибли- жением решения.
Расчёт выполняем по формуле:
X
i
= X
i-1
+ Y
i-1
/dY
i-1
, где индекс i – указы- вает на номер строки в таблице решения задачи. Однако, как видно из рис. 12 надо учитывать направление дви- жения к следующему приближе- нию (по направлению оси «Х» или против него, как показывают стрелки). Для определения направ- ления можно воспользоваться зна- чениями функции и её производ- ной. Знак при их произведении по- казывает направление движения.
Если произведение этих чисел от- рицательно, то двигаемся право по оси и наоборот при положитель- ном произведении двигаемся влево от больших значений Х к меньшим. Данные подход будет использован и при реализации ко- нечно-разностных методов для определения знака шага поиска (h).
Рис. 11. Определение следующего приближения для поиска
Y=F(X)
Y'=F(X)
Y
0
Y
0
X
0
X
1
Y
X
0
X
0
X
0
Y>0, dY<0
Y<0, dY>0
Y>0, dY>0
Y<0, dY<0
H
-H
Рис. 12. Определение направле- ния поиска следующего прибли- жения

18
В результате можно записать общую формулу для вычисления сле- дующего приближения во второй строке таблицы для столбца Х: если произведения значений функции на её производную является отрица- тельным числом, то смещение (Y
i-1
/dY
i-1
) прибавляется иначе отнима- ется. В таблице эта формула выглядит так:
=D8+ЕСЛИ(E8*F8<0;
ABS(E8/F8); -ABS(E8/F8))
. Функция Abs() убирает знак этого произведения, а основной знак стоит перед данной функцией.
Копируем две имеющиеся формулы (расчёт функции и её производ- ной) из первой строки таблицы во вторую строку. Теперь вторую строку таблицы можно копировать вниз, пока значения функции не прибли- зятся к нулю (рис. 13а). Для контроля за сходимостью решения можно достроить таблицу столбцом «
Erf
» (ошибка решения) с формулой:
=Abs(D8-D9)
, где
D8
– значение функции из предыдущей строки;
D9
– зна- чение функции из текущей строки. Результаты показаны на рис. 13б. а. б.
Рис. 13. Результаты расчёта: а. – без контроля; б. – с контролем приближения
Построим траекторию поиска при решении этой задачи. Готовим таблицу вводя её заголовок и шапку для столб- цов (рис. 14) и начинаем её заполнять данными из расчётной таблицы. В первую строку ставим ссылки на начальную точку поиска (X, Y из пер- вой строки таблицы с решением задачи). Из этой точки по касательной (прямая линия) надо пере- меститься в следующее значение Х, значит в столбец «
Х
» ставим ссылку на второе приближение, а в столбец «
Y
» зна- чение 0 (пересечение с осью «X»). В третей строке повторяем второе приближение по «Х», а в столбец «
Y
» ставим расчётное значение функ- ции для него из таблицы решения задачи.
Рис. 14. Таблица для построения траектории

19
Далее продолжаем по той же схеме. В резуль- тате на четыре шага поиск имеем следующую таблицу (рис. 15).
Добавляем эти данные в диаграмму командой
«
Выбрать данные
» – «
Добавить
» и получаем следу- ющий график (рис. 16). Мы получили лист, кото- рый позволяет реализовывать поиск корня урав- нения и графически наблюдать за решением за- дачи. Можно изменять начальное приближение и график будет перестраиваться автоматически.
Как видим, начальное приближение не всегда приводит к решению задачи, процедура поиска может и расходиться, если за- дать начальное приближение в близи минимума функции
(7,5). Для других функций вы- бор начального приближения может быть ещё более жёст- ким.
Попробуем реализовать по- иск корня для уравнения с не- известной производной, вос- пользовавшись методом хорд.
Для построения хорды не бу- дем задавать вторую точку, а вычислим дополнительную точку на расстоянии 0,1*Х от заданной вперёд по оси Х.
Так как решение поиска корня подобны, поэтому сделаем копию со- зданного листа с полиномом в качестве уравнения, заменим заголовок на «
Произвольная функция
» и уберём строки со значениями коэффициен- тов полинома. В качестве исследуемой функции используем
Y=30-0,8 X
2
+5Ln(X
3
) и запишем рядом с заголовок листа (рис. 17).
Рис. 17. Вид заголовка листа.
Рис. 15. Таблица траектории поиска
Рис. 16. График функции с траекторией поиска корня
-10
-5
0
5
10
15
20
25
0
2
4
6
8
10
Y=F(X)
Траектория поиска

20
Исправляем уравнение в столбце
Y=F(X)
таблице «Табулирование функции» на следующую фор- мулу:
=30-0,8*A5^2+5*LN(A5^3)
, где:
A5
- ссылка на значение Х. По- том копируем формулу на весь столбец. В первой строке таблицы возникает ошибка
#ЧИСЛО!
, свя- занная с отсутствием решения в точке «0», исправляем значение параметра Х на 0,1. График стро- ится автоматически и имеет вид
(рис. 18).
В таблицу
«
Решение задачи
»
в первую ячейку столбца «
Y=F(X)
» копируем введённую формулу функции из первой таблицы, стол- бец производной нам не нужен, удаляем весь столбец, выделив его с имени заголовка столбца. Во вто- рую строку для приближения за- писываем формулу для расчёта следующего приближения. Разбе- рём эту формулу более подробно, вычислять следующее приближе- ние будем через подобие треуголь- ников (рис. 19). От заданной точки
(X, Y) строим новую точку с шагом
dX. Имеем два треугольника: пер- вый вокруг двух точек и второй от заданной точки до оси Х, они по- добны и катеты первого треугольника известны, а во втором один из катетов надо определить. Получаем формулу:
i
dX
i
i
i
i
i
i
i
dX
Y
Y
dX
Y
X
X
Y
X
X
Y
Y
dX











1 1
;
(1)
Как можно выбрать шаг по Х? Можно его задать конкретным чис- лом, но решение в этом случае не будет сходиться из-за большого шага
Рис. 18. График произвольной функции
-20
-10
0
10
20
30
40
50
0
2
4
6
8
10
Y=F(X)
Рис. 19. Расчёт приближения в методе хорд
Y
i
Y
+dX
X
i
X
+dX
dX
(Y
i
- Y
+dX
)
X
i+1
- X
i
X
i+1

21 для второй точки вблизи корня. Шаг должен быть пропорционален зна- чению Х, например, 0,1*Х, такой подход имеет одну специфическую точку Х=0, где надо принять шаг равный конкретному числу, например,
0,1. Для дополнительной точки использовать множитель 1,1 нельзя, если возможны отрицательные значения Х, поэтому координата Х до- полнительной точки вычисляются по формуле: Х
i
+ Abs(X
i
)*0,1.
Окончательный вид формулы показан далее:
=D5-ЕСЛИ(D5<>0; ABS(D5)*0,1*(30-0,8*(D5)^2+5*LN((D5)^3))/((30-
0,8*(D5+ABS(D5)*0,1)^2+5*LN((D5+ABS(D5)*0,1)^3))-(30-
0,8*(D5)^2+5*LN((D5)^3)));0,1*(30-0,8*(D5)^2+5*LN((D5)^3))/((30-
0,8*(0,1)^2+5*LN((0,1)^3))-(30-0,8*(D5)^2+5*LN((D5)^3))))
Разберём её более подробно. В операторе «
Если(
» проверяем входной параметр
D5
на его неравенство нулю, в этом случае строим формулу:
(Y(X)*dX)/(Y(X+dX)-Y(X)), где X+dX – заменяем выражением Х
i
+
Abs(X
i
)*0,1. В случае равенства входного параметра «0» просто берём вторую точку с координатой 0,1 по Х. Полученную формулу растяги- ваем до конца таблицы и видим решение задачи и график с траекторией поиска (рис. 19).
Рис. 19. Расчётная таблица и график
Можно реализовать решение и через вычисление производной для заданной функции через конечно-разностный её аналог, что приведёт к подобной сложной формуле:
=ЕСЛИ(D5<>0;((30-0,8*(D5+ABS(D5)*0,01)^2+5*LN((D5+ABS(D5)*0,01)^3))-
(30-0,8*(D5-ABS(D5)*0,01)^2+5*LN((D5-
ABS(D5)*0,01)^3)))/(ABS(D5)*0,02);(30-0,8*(0,01)^2+5*LN((0,01)^3))-(30-
0,8*(-0,01)^2+5*LN((-0,01)^3))/(0,02))
-40
-30
-20
-10
0
10
20
30
40
50
0
2
4
6
8
10
12
Траектория поиска
Y=F(X)

22
Где производная вычисляется по центральной разности с помощью двух точек, равноудалённых на ±0,01*Х от заданной точка.
X
Y
Y
X
Y
dX
dY
X
X











2
(2)
Как и в предыдущей формуле имеется одна исключаемая точка при
Х=0, где шаг задаём жёсткой равным ±0,01. В результате получаем тот же результат расчёта, что показан на рис. 19.
Для упрощения этих задач целесообразно воспользоваться функци- ями пользователя, созданными в среде языка Бейсик (VBA for Excel).
Требуется создать две функции: саму функцию Y=F(X) для которой ищем корень и функцию вычисления производной для заданной функ- ции Y=F(X).
Реализуем данного ре- шение. Вызываем редак- тор
VBA командой
Alt+F11
, открывается окно редактора (рис. 20).
Командой «
Insert
» –
«
Module
» создаём модуль программы (рис. 21).
Рис. 21. Окно с созданным модулем
Находясь в его (модуля) поле создаём функцию командой «
Insert
» –
«
Procedure
», которая открывает окно (Рис. 22) и в нем записываем имя функции, например, «
My_fun01
» и выбираем тип процедуры как
«
Function
». Завершаем ввод информации кнопкой «
Ok
». На панели мо- дуля появляется заготовка процедуры,
Public Function My_fun01()
End Function
Рис. 20. Окно редактора VBA