Численное дифференцирование

Раздел № 5

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

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

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

Пусть известен график функции у = f (х ) на отрезке [а ,b ].Можно построить график производной функции, вспомнив ее геометрический смысл. Воспользуемся тем фактом, что производная функции в точке х равна тан­генсу угла наклона к оси абсцисс касательной к ее графи­ку в этой точке.

Если х = х 0 ,найдем у 0 = f (x 0)с помощью графика и затем проведем касательную АВ к графику функции в точке (х 0 , y 0) (рис. 5.1). Проведем прямую, параллельную касательной АВ, через точку (-1, 0) и найдем точку у 1 ее пересечения с осью ординат. Тогда значение у 1 равно тан­генсу угла наклона касательной к оси абсцисс, т. е. про­изводной функции f (x )в точке х 0:

у 1 = = tgα = f ¢ (x 0), и точка М 0 (х 0 , у 1) принадлежит графику производной.

Чтобы построить график производной, необходимо разбить отрезок [а , b ]на несколько частей точками х i , затем для каждой точки графически построить значение производной и соединить полученные точки плавной кри­вой с помощью лекал.

На рис. 5.2 показано построение пяти точек М 1, М 2 ,... , М 5 и графика производной.

Алгоритм построения графика производной:

1. Строим касательную к графику функции у = f (x )в точке (х 1 , f (x 1));из точки (-1, 0) параллельно касатель­ной в точке (х 1 , f (x 1)) проведем прямую до пересечения с осью ординат; эта точка пересечения дает значение про­изводной f ¢ (х 1).Строим точку М 1 (х 1 , f ¢ (х 1)).

2. Аналогично построим остальные точки М 2 , М 3 , М 4 и М 5 .

3. Соединяем точки М 1 , М 2 , М 3 , М 4 , М 5 плавной кривой.

M 4

Полученная кривая является графиком производной.

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

Замечание . Если в алгоритме построения графика производ­ной вместо точки (-1, 0) взять точку (-l ,0), где l > 0, то график будет построен в другом масштабе по оси ординат.

5 . 2 .Разностные формулы

а) Разностные формулы для обыкновенных производных

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

y i = f (x i ), x i = a+ ih , i = 0, 1, ... , n; h =

Мы рассматриваем случай равномерного распределения точек на отрезке [a , b ]. Для приближенного вычисления производных в точках x i можно использовать следующие разностные формулы , или разностные производные .

Так как предел отношения (5.1) при h ® 0 равен пра­вой производной в точке х i , то это отношение иногда на­зывают правой разностной производной в точке x i .По аналогичной причине отношение (5.2) называют левой разностной производной в точке x i .Отношение (5.3) на­зывают центральной разностной производной в точке x i .

Оценим погрешность разностных формул (5.1)–(5.3), предполагая, что функция f (x ) разлагается в ряд Тейло­ра в окрестности точки x i :

f (x ) = f (x i )+ . (5.4)

Полагая в (5.4) х = x i + h или х = х i - h , получим

Непосредственной подстановкой разложений (5.5) и (5.6) в формулу (5.10) можно получить зависимость между второй производной функции и разностной формулой для производной второго порядка .

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

На рис. 5.8, а изображена кривая полученная экспериментально на установке (рис. 5.6). Определение углового ускорения (искомой функции) проводят графическим дифференцированием по соотношению:

(5.19)

Тангенс угла наклона касательной к кривой в некоторой точке i представляют в виде отношения отрезков , где К – выбранный отрезок интегрирования (рис. 5.8, б )

После подстановки этого соотношения в соотношение (5.19) полу­чают

где - ордината искового графика углового ускорения;

Масштаб искомого графика ; единицы СИ: = мм; = мм/(рад с -2).

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

Графическое дифференцирование рассмотренным методом каса­тельных имеет относительно низкую точность. Более высокую точность получают при графическом диф­ференцировании методом хорд (рис. 5.8, в и г ).

На заданной кривой отмечают ряд точек 1 ", 2 ", 3" , которые соединяют хордами, т.е. заменяют заданную кривую ломаной ли­нией. Принимают следующее, допущение: угол наклона касательных в точках, расположенных посередине каждого участка кривой, ра­вен углу наклона соответствующей хорды. Это допущение вносит некоторую погрешность, но она относится только к данной точке. Эти погрешности не суммируются, что обеспечивает приемлемую точность метода.

Остальные построения аналогичны ранее описанным при графи­ческом дифференцировании методом касательных. Выбирают отре­зок (мм); проводят лучи, наклоненные под углами до пересечения с осью ординат в точках 1 ", 2 ", 3 " ... , которые переносят на ординаты, проведенные в середине каждого из интервалов. Полученные точки 1 *, 2 *, 3 * являются точками иско­мой функции .

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

Дифференцирование функции f(x) , заданной (либо вычисленной) в виде массива чисел, выполняют методом численного дифферен­цирования с применением ЭВМ.

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

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

где Δх – очень малая конечная величина.

При достаточно малых значениях Δх, можно с приемлемой точностью получить величину производной функции в точке. Для вычисления производной в MS Excel будем использовать приведенную выше формулу. Рассмотрим технологию вычисления производной на примере .

Пример 1.18 Найти производную функции у = 2х 3 + х 2 в точке х=3. Заметим, что производная приведенной функции в точке х=3, вычисленная аналитическим методом, равна 60 - это значение нам понадобится для проверки результата, полученного путем вычисления численным методом.

Задачу вычисления производной в табличном процессоре можно решать двумя способами.

Решение первым способом

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

2*А2 ^ 3+А2 ^ 2.

Зададим окрестность точки х=3 достаточно малого размера, например значение слева х k =2,9999999, а значение справа х k +1 =3,00000001, и введем эти значения в ячейку А2 и А3 соответственно. В ячейку С2 введем формулу вычисления производной =(В3-В2)/(А3-А2).

В результате вычисления в ячейку С2 будет выведено приближенное значение производной заданной функции в точке х=3, величина которой равна 60, что соответствует результату, полученному аналитически (рис.1.24).

Решение вторым способом

Введем в ячейку рабочего листа А2 заданное значение аргумента, равное 3, в ячейке В2 укажем достаточно малое приращение аргумента - (1E - 9), в ячейку С2 введем формулу для вычисления производной

=(2*(А2+В2) ^ 3+(А2+В2) ^ 2-(2*А2 ^ 3+А2 ^ 2))/В2.

После нажатия клавиши получим результат вычисления 60,0000.

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

Вычисление локальных экстремумов функции

Напомним, что функция Y=f(x) имеет экстремум при значении х = х k если производная функции в этой точке равна нулю.

Если функция f(x) непрерывна на отрезке [а, b] и имеет внутри этого отрезка локальный экстремум, то его можно найти, используя надстройку Excel Поиск решения.

Рассмотрим последовательность нахождения экстремума функции на примере

Пример 1.19 Задана неразрывная функция у = х 2 + х + 2. Требуется найти ее экстремум (минимальное значение) на отрезке [-2; 2].

Решение

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

В ячейку В3 введем формулу, определяющую заданную функциональную зависимость. Вместо переменной х в этой формуле должна быть ссылка на ячейку А3: =А3^2+A3+2.

Выполним команду меню Сервис/Поиск решения.

В открывшемся окне диалога Поиск решения в поле Установить целевую ячейку укажем адрес ячейки, содержащей формулу (В3), установим переключатель Минимальному значению, в поле Измени ячейки укажем адрес ячейки, в которой содержится переменная х-A3.

Добавим два ограничения в соответствующее поле: A3 > = - 2 и A3<=2 (рис. 1.25).

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

Щелкнем на кнопке Выполнить. В ячейке А3 будет вычислено значение аргумента х функции, при котором она принимает минимальное значение, а в ячейке В3 – минимальное значение функции.

В результате выполнения вычислений в ячейке А3 будет получено значение независимой переменной, при котором функция принимает наименьшее значение -0,5, а в ячейке В3 – минимальное значение, равное 1,75.

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

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

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

Используя приведенную технологию численного вычисления производной функции в заданной точке, проверим, является ли найденная точка х = -0,5 точкой экстремума функции у = х 2 + х + 2. Решение приведено на рисунке.

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

Пример 1.20 Требуется найти значения аргумента в диапазоне [-1; 1], при которых функция у = х 2 + х + 2 имеет экстремумы.

Решение

Табулируем заданную функцию с шагом 0,2.

Применяя второй из приведенных способов вычисления производной, вычислим значения функции у = f(x + dx).

Вычислим значения производной при каждом табличном значении аргумента.

Анализируя полученные значения производных функции в точках, находим, что производная меняет знак в интервале значений аргумента (-0,6;-0,4), следовательно, на этом интервале есть точка экстремума. Кроме того, заметим, что знак производной меняется с минуса на плюс, следовательно, точка экстремума является минимумом функции.

Применяя инструмент Подбор параметра или Поиск решения для решения уравнения Y(x) = 0

относительно х, вычислим точное значение аргумента, при котором исходная функция принимает экстре малыше значение (-0,5) (рис. 1.26).

Полученное значение производной исследуемой функции в точке х =-0,5 равно нулю, следовательно, в этой точке функция имеет экстремум.

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

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

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

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

Во второй части статьи на «живом» примере рассмотрим вычисление производной по приближенной формуле численного дифференцирования с применением выражений в конечных разностях и разберемся в вопросе – можно ли используя приближения производных конечными разностями по прогибам балки определять действующие в сечениях нагрузки?

Минимум теории.

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

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

y’ (x )=lim (Δy /Δx ) при Δx →0

Геометрический смысл производной функции в точке – это тангенс угла наклона к оси x касательной к графику функции в этой точке.

tg (α )=Δy /Δx

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

y’ (x ) i ≈(Δy /Δx) i =(y i +1 -y i -1 )/(x i +1 -x i -1 )

Конечными разности называют потому, что они имеют конкретное, измеримое, конечное значение в отличие от величин, стремящихся к нулю или бесконечности.

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

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

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

Дано:

На балку длиной 8 метров с шарнирными опорами по краям изготовленную из двух спаренных стальных (Ст3) двутавров 30М опираются 7 прогонов с шагом 1 метр. К центральной части балки крепится площадка с оборудованием. Предположительно усилие от покрытия, передаваемое через прогоны на балку, во всех точках одинаково и равно F 1 . Подвесная площадка имеет вес 2*F 2 и крепится к балке в двух точках.

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

На рисунке ниже показана расчетная схема задачи и общий вид эпюр.

На следующем скриншоте представлены исходные данные.

Расчетные исходные данные:

3. Погонная масса двутавра 30М:

γ =50,2 кг/м

Сечение балки составлено из двух двутавров:

n =2

Удельный вес балки:

q =γ *n *g =50,2*2*9,81/1000=0,985 Н/мм

5. Момент инерции сечения двутавра 30М:

I x1 =95 000 000 мм 4

Момент инерции составного сечения балки:

I x =I x 1 *n =95 000 000*2=190 000 000 мм 4

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

R =(q *z max +8*F 1 +2*F 2 )/2=(0,985*8000+8*9000+2*50000)/2=85 440 Н

В расчете учитывается собственный вес балки!

Задача:

Найти значения изгибающего момента M xi в сечениях балки аналитически по формулам сопротивления материалов и методом численного дифференцирования расчетной линии прогибов. Сравнить и проанализировать полученные результаты.

Решение:

Первое, что мы сделаем, это выполним расчет в Excel поперечных сил Q y , изгибающих моментов M x , углов поворота U x оси балки и прогибов V x по классическим формулам сопромата во всех сечениях с шагом h . (Хотя, в принципе, значения сил и углов нам в дальнейшем не понадобятся.)

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

Использованные в расчетах формулы можно посмотреть .

Итак, нам известны точные значения моментов и прогибов.

Из теории мы знаем, что:

Угол поворота – это первая производная прогиба U =V’ .

Момент – это вторая производная прогиба M =V’’ .

Сила – это третья производная прогиба Q =V’’’ .

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

M xi =V y ’’ ≈((V i +1 -2*V i +V i -1 )/h 2)*E *I x

Итог расчетов мы видим в ячейках M5-M54.

Точные значения моментов, рассчитанные по аналитическим формулам сопромата с учетом веса самой балки, отличаются от найденных по приближенным формулам вычисления производных незначительно. Моменты определены весьма точно, судя по относительным погрешностям, рассчитанным в процентах в ячейках N5-N54.

ε =(M x -V y ’’ )/M x *100%

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

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

Однако...

Увы, не стоит думать, что на практике легко получить необходимые высокоточные результаты измерений прогибов сложно нагруженных балок!

Дело в том, что измерения прогибов требуется выполнять с точностью ~1 мкм и стараться максимально уменьшать шаг замеров h , «устремляя его к нулю», хотя и это может не помочь избежать ошибок.

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

Сегодня есть приборы — лазерные интерферометры, обеспечивающие высокую скорость, стабильность и точность измерений до 1 мкм, программно отсеивающие шум, и еще много чего программно умеющие, но их цена – более 300 000$...

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

Если раньше максимальная ошибка не превышала 0,7%, то сейчас (в сечении i =4) превышает 23%, хотя и остается приемлемой в наиболее опасном сечении (ε 21 =1,813%).

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

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

Вычисление производной численными методами всегда необходимо выполнять очень осторожно!

Уважаемые читатели, отзывы и комментарии к статье, размещайте в специальном блоке ниже статьи.

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

Прошу УВАЖАЮЩИХ труд автора скачать файл с примером ПОСЛЕ ПОДПИСКИ на анонсы статей.

Похожие статьи