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

23 которую надо наполнить содержанием. Сначала в список формальных параметров функции, расположенный внутри скобок после её имени помещаем единственные параметра «
Х
».
Потом в тело функции (пустая строка между заголовком и концом проце- дуры) записываем необходимый код про- граммы. В нашем случае приравниваем имени функции расчёт нашей функции с использование указанного параметра:
Public Function My_fun01(X)
My_fun = 30 - 0.8 * X ^ 2 + 5 * Log(X ^ 3)
End Function
Вторая функция, которая нам по- надобится в дальнейшем – это вычисление производной для любой функции. Вычис- ление производной будем выполнять чис- ленным методом. Существует достаточно много формул для вычисления производ- ной в заданной точке исследуемой функ- ции. Рассмотрим самый простой из них
(рис.1.14). От заданной точки Х
0
отступим на достаточно малую величину Δ по Х вправо и влево и найдём в этих точках зна- чения функции и вычислим производную по формуле (2).
Надо помнить, что отнимать надо все- гда от правой точки левую. Чем меньше шаг, тем выше точность реше- ния. Как можно опередить шаг? Если мы ведём вычисления с точностью до 5%, то шаг должен быть не более 2% от значения аргумента Х. При- нимаем шаг равный 1% от значения параметра, что приводит к следую- щему алгоритму расчёта:








X
X
X
X
01
,
0
;
0 01
,
0
;
0
(3)
Если Х=0, принимаем шаг равным 0,01, в других случаях шаг равня- ется 0,01 от значения Х.
Рис. 22. Окно создания процедуры
Y
X
0
X
0
Y
0
X
-

X
+

Y
-

Y
-

1
2
3
Рис.1.14. Вычисление производной

24
Теперь создадим эту функцию, назвав её «
Difr_Ur
». Для этого встав- ляем новую процедуру с этим именем и вводим следующий код:
Public Function Difr_Ur(Xt, Fn As String)
If Xt <> 0 Then dH = 0.01 * Xt
Else dH = 0.01
End If
XL = Application.Run(Fn, X - H)
XR = Application.Run(Fn, X + H)
Difr_Ur = (XR - XL) / (2 * dH)
End Function
Функции пользователя получены и их можно использовать на рабо- чем листе Excel. Копируем лист со столбцом для вычисления производ- ной и изменяем формулу для расчёта «
Y=F(X)
» в первой и второй табли- цах с использованием созданной функции пользователя:
=My_Fun01(D5)
Для вызова данной функции открываем мастер построения функций и выбрать категорию «
Определённые пользователем
», там находим нашу функцию или просто набрать «
=My_
» и она появится в списке возмож- ных функций. Указываем ссылку на фактический параметр функции
(ячейка в столбце Х) и завершаем ввод кнопкой «
Ok
». Копируем фор- мулу до конца таблицы.
Для вычисления производной воспользуемся второй созданной функцией «
Difr_Ur
», которую запишем в соответствующий столбец вто- рой таблицы. Данная функция имеет два формальных параметра
Diff_Ur(D5;"My_Fun01")
:
D5
– значение Х и "My_Fun01"
– имя функции, про- изводную от которой хотим найти. Второй параметр записывается как текстовая строка в самой функции или ссылка на ячейку, где заранее можно вбить наименование функции. Копируем формулу до конца таб- лицы и расчёт работает.
Теперь мы имеем рабочий лист, на котором достаточно сменить имя новой функции и расчёт будет выполняться автоматически. Попробуем исследовать функцию, у которой легко найти области существования решений данным методом. В качестве рабочей функции воспользуемся обратным тангенсом (arctan). Создадим новую функцию пользователя
«
My_fun1
» со следующей формулой:
10 * Arctg(X - 5) – 8
, график которой показан на рис. 23. Такая кривая при использовании этого метода очень чувствительна к выбору начальной точки, проверим это.

25
Корень функции лежит около 6, зададим начальную точку равную 7,7 и получим решение (рис. 24а), а при начальном приближении 7,8 решения нет, и процедура поиска расходится
(рис. 24б), что так же видно из рас- чётной таблицы по столбце «Erf» в котором в первом случае значения стремятся к нулю, а во втором к боль- шим числам. а. б.
Рис. 24. Поиск корня методом касательных: а. – решение сходится; б. – решение расходится
Второй проекционный метод, который мы рассмотрим – это метод итераций.
-10
-5
0
5
10
15
20
25
0
2
4
6
8
10
12
Y=F(X)
Траектория поиска
-10
-5
0
5
10
15
20
25
-10 -8 -6 -4 -2 0
2
4
6
8 10 12
Y=F(X)
Траектория поиска
Рис. 23. Функция арктангенса
-10
-5
0
5
10
15
20
25
0
2
4
6
8
10
12
Y=F(X)

26
0>0>

27 и произвольное уравнение. Создаём новый лист Excel и формируем в верхней его части данные с коэффициентами полинома (см. рис. 6). Да- лее строим таблицу нашей функции, чтобы видеть её график, как делали в прошлой задаче. Для наблюдения за ходом решения надо построить ещё две функции X
0
=F(X) и X
1
=F(X), ко- торые добавляем в эту же таблицу (рис.
26) и для столбца «
X0=F(X)
» записываем формулу =Х (
=A8
), а для столбца
«
X1=F(X)
» используем преобразованную функцию от полинома, выразив линей- ный коэффициент уравнения, в нашем случае уравнение имеет вид:
=-($B$2+$B$4*C8^2)/$B$3
, где
$B$2
,
$B$3
,
$B$4
– коэффициенты полинома;
C8
– ар- гумент Х. Копируем формулы до конца таблицы и строим два графика: Y=F(X) (рис. 27а), X
0
=F(X) и Х
1
=F(X)
(рис. 27б). а. б.
Рис. 27. Графики исследуемой функции: а. – сама функция от аргумента; б. – начальный и следующий Х от аргумента.
В данном случае поиск ведётся от линии Х
0
=F(X), которая является диагональю на графике (рис. 27б).
Готовим таблицу для решения задачи, она содержит три основных столбца:
- столбец «
Х0
» с заданным в первой строке начальным приближением и полученным следующим значением из предыдущей строки;
- столбец «
Х1
» с функцией для вычисления следующего приближения от начального;
-10
-5
0
5
10
15
20
25
0
2
4
6
8
10
Y=F(X)
0
2
4
6
8
10
0
2
4
6
8
10
X0=F(X)
X1=F(X)
Рис. 26. Фрагмент таблицы

28
- столбец «
Erf
» с разностью приближений (Х
0
-Х
1
) для контроля реше- ния задачи.
Формулы в таблице «Решение задачи» и её окончательный вид по- казаны на рис. 28. а. б.
Рис. 28. Таблица «Решение задачи»: а. – формулы; б. – результат
Таблица «Траектория движения» строится от начального приближения в первой строке (Х
0
,
Х
0
), потом к точке (Х
0
, Х
1
) во второй строке. В третьей строке повторяется первая с ссылками на второе приближение и т.д. Данная таблица имеет вид (рис. 29). Сравните данные из таблиц
(рис. 28б) и (рис. 29), чтобы понять принцип по- строения.
По данным этой таблицы добавляем на вто- рой график (рис. 27б) траекторию поиска ответа, как это делалось раньше (рис. 30).
Реализуем данное решение для трансцендентного уравнения вида:
Y = 0,005*X^3+Ln(0,1*X)
Одним из вариантов преобразо- вания уравнения можно рассматри- вать следующий вид:
X1=(-LN(0,1*X0)/ 0,005)^(1/3)
Скопируем лист с методом ите- раций и заменим в нем функции для расчёта Y и Х
1
и получим реше- ние задачи, где показан поиск ре- шения на графике функции и поиск корня по траектории поиска (рис.
31б).
Рис. 29. Таблица траектории поиска
Рис. 30. Поиск корня.
0
2
4
6
8
10
0
2
4
6
8
10
Траектория поиска
X0=F(X)
X1=F(X)

29 а. б.
Рис. 31. Графики поиска корня функции
В данном методе использовать программы в VBA не имеет смысла, из-за простоты записи функций и их индивидуальности, что исключает возможность упрощения данного решения через создание макросов или функций пользователя.
Для реализации конечно-разностной схемы решения разберём реше- ние задачи нахождения корня методом деления отрезка пополам.
Метод деления отрезка пополам (дихотомии)
Задаётся интервал, в котором имеется корень (проверка наличия корня внутри интервала выполняется вычислением значений функции на его границах с последующей проверкой этих значений на разные знаки через операцию умножения, если произведение меньше нуль, то знаки на границах интервала разные). Затем интервал разбивается на два, и выбирается тот, в котором остался корень. Процедура повторя- ется, пока интервал не станет меньше заданной погрешности erf.
Алгоритм метода показан на рис. 32. Обычно проверяется левый ин- тервал, а после проверки выбирается нужный. Ответом является сере- дина последнего интервала. Для ускорения работы алгоритма можно проверять и попадание в сам корень, в этом случае произведение стано- вится равным нулю, а корень оказывается на одной из границ интер- вала. Разберём данную процедуру подробнее:
 задаём уравнение Y=f(X), погрешность вычисления erf и левую гра- ницу Х
л
, проверяем наличие корня на левой границе интервала. Если он есть, завершаем работу.
-15
-10
-5
0
5
10
0
2
4
6
8
10
Поиск
Y=F(X)
0
2
4
6
8
10
12
14
0
2
4
6
8
10
Траектория поиска
X0=F(X)
X1=F(X)

30
 задаём правую границу Х
пр
, проверяем наличие корня на правой границе интервала, если он есть, завершаем работу про- граммы. Проверяем произведе- ния f(X
л
)·f(X
пр
).Если значение больше нуля, то корня в интер- вале нет. Повторяем ввод пра- вой границы. Если произведе- ние меньше нуля, то переходим к поиску корня.
 находим середину интервала
Х
ср
=(Х
пр
+Х
л
)/2, проверяем наличие корня в данной точке, если он есть, завершаем работу программы. Проверяем значе- ние произведения для левой по- ловины интервала f(X
л
)·f(X
ср
).
Если значение меньше нуля, то левая половина содержит ко- рень и правую часть интервала можно отбросить, приравняв
Х
пр
= X
ср
, иначе корень нахо- дится в правой части, и тогда отбрасываем левую Х
л
= X
ср
 процедура повторяется до тех пор, пока (Х
пр
-Х
л
)>erf. Ответом является среднее значение по- следнего найденного интер- вала.
Реализуем это решение на ли- сте Excel. Создаём новый лист и пишем заголовок метода «Метод деления отрезка пополам».
Ниже готовим шапку таблицы с тремя сдвоенными столбцами
Рис. 32. Алгоритм метода деле- ния отрезка пополам (дихото- мии)
Метод
деления отрезка
пополам
(
дихотомии)
Исходные
данные
Y=f(X) X
л
,
erf
Xcp=(X
пр+Xл)/2
Abs(X
пр
-X
п
)>erf
Ответ
(X
пр+Xл)/2
Конец
X
пр
f(X
л)·f(Xпр)>0
Да
Нет
f(X
л)·f(Xcр)>0
X
л
=X
ср
X
пр
=X
ср
Да
Нет

31 для координат левой границы интервала, его середины и правой границы интервала. В следующей строку указываем имена коор- динат «Х» и «Y». Завершаем шапку таблицы столбцом «Ошибка».
Окончательный вид заголовка и шапки показан на рис. 33.
Рис. 33. Заголовок метода и шапка таблицы.
При оформлении шапки таблицы объединяем по две ячейки по горизонтали для границ интервала и середины, для столбца оши- бок объединяем две ячейки по вертикали.
Теперь можно начинать заполнять таблицу данными, сначала отмечаем координаты «Х» левой и правой границ фоном для напоминания о необходимости ввода данных в эти ячейки. Потом заполняем все столбцы «Y» нужной функцией и убеждаемся в наличии корня в заданном интервале (значения функций на левой и правой границах должны иметь разные знаки). Теперь вводим формулу расчёта середины интервала =(X
R
+ X
L
)/2. В столбец ошибки вносим формулу =Abs(X
R
– X
L
). Первая строка заполнена формулами, переходим ко второй строке.
Сначала надо заполнить формулы для столбцов «Х» на грани- цах интервала через условие ЕСЛИ() и наличия разных знаков у функций (Y) для левой границы и середины. Проверку лучше ве- сти по одному и тому же условия, чтобы не попасть в разночтения условий для каждого интервала. Формула должна звучать следу- ющим образом: если Y
L
и Y
C
имеют разные знаки, то левая граница сохраняет своё значение, а правая меняет своё на среднее значе- ние и наоборот, при одинаковых знаках Y
L
и Y
C
X
L
становится рав- ным X
C
, а X
R
сохраняет своё значение. Отсюда можно построить две формулы для X
L
и X
R
:
=ЕСЛИ(B4*D4<=0;A4;C4) =ЕСЛИ(B4*D4<=0;C4;E4)
где:
B4 и
D4
– значения Y на левой границе интервала и в середине;
A4
,
C4
,
E4
– значения координат Х на левой границе, в середине и на правой границе соответственно.