Модели и методы решения нелинейных уравнений

Моделирование динамических систем: решение нелинейных уравнений

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

На выходе мы получили траекторию полета снаряда, что дает нам ориентировочные представления о характеристиках орудия: при заданных параметрах мы получили дальность стрельбы чуть более 2,5 км и высоту подъема снаряда чуть выше 800 метров. Точнее мы сказать не можем, вернее можем, если с карандашиком по сетке определим координаты нужных точек на графике. Но это, как известно, не наш метод. Хорошо бы получить эти данные с точностью, обеспечиваемой используемыми нами инструментами и без ручного труда. Вот об этом мы сегодня и поговорим.

1. Постановка задачи

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

как функции времени. Высота полета снаряда это y(t). Если мы определим в какой момент времени высота становится равна нулю, то есть решим уравнение

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

Как найти максимальную высоту подъема снаряда? Из графика траектории видно, что по мере подъема снаряда его траектория становится всё более и более пологой, пока в экстремальной точке скорость на мгновение становится горизонтальной и дальше снаряд движется вниз.

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

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

Просто? Более чем.

2. Уравнение, которого нет

И вот тут котенка Гава, как известно, ждут неприятности. Начнем с того, что даже если уравнение задано в виде формулы (аналитически) не всегда удается найти его решение. Вот например

Как вам? Простенько, но попробуйте найти икс, используя всё то, чему вас учили в школе. То-то же…

Введем замену , тогда
\begin
&u \, e^ <2 — u>= 1 \\
&u \, e^ <-u>= e^ <-2>\\
&-u \, e^ <-u>= -e^ <-2>
\end
Пусть теперь , тогда

Теперь делаем финт ушами. Математики прошлого хорошо поработали за нас. Если задана функция вида

то обратная ей функция, называется W-функцией Ламберта.

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

откуда, раскручивая назад все замены получаем ответ

Приближенно этот ужас равен

Решение такого уравнения придется искать численно, тем более что очевидно его графическое решение

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

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

Создадим новый скрипт в том же каталоге, где расположены файлы ballistics.m, f.m и f_air.m. Назовем его, например cannon.m. Для начала зададимся параметрами снаряда, начальной скоростью и направлением выстрела

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

Обратите внимание, теперь в качестве моментов времени, передаваемых функции решения уравнения мы берем всего два значения: начальный момент времени (t = 0) и интересующий нас момент времени. Соответственно, переменная solv будет содержать два вектора фазовых координат: начальный и тот который нам нужен. Собираем все компоненты конечной точки фазовой траектории в вектор Y и возвращаем его значение из функции.

Теперь нам ничего не стоит определить зависимость высоты полета снаряда от времени

протестируем полученную функцию

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

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

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

3. Принципы численного решения нелинейных уравнений

Методы численного решения нелинейных и трансцендентных уравнений заточены под решение уравнений вида

К такой форме легко привести любое уравнение. Например

где . Корни этого, эквивалентного уравнения, равны корням исходного. Если мы построим график функции f(x), то увидим такую картинку

корни уравнения это значения аргумента в тех точках, где график пересекает ось x.

Все методы численного решения таких уравнений включают в себя два этапа:

1. Поиск начального приближения корня
2. Уточнение корня с заданной погрешностью

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

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

и выполняю расчеты по рекуррентной формуле

В нашем примере

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

Ага, видно, что начиная с шестой итерации четвертый знак получающегося числа остается неизменным. Значит мы нашли корень уравнения с погрешностью менее 10 -4 на шестой итерации.

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

Ограничимся членами первого порядка малости, заменив саму функцию f(x) касательной к её графику в точке x₀

и приравняв полученное выражение к нулю, решим его относительно x

Получаем итерационную формулу

В нашем примере

В качестве начального приближения берем и пытаемся выполнять итерации

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

Эти примеры я привел, чтобы объяснить общий принцип. Кроме этих двух методов существует ещё масса методов, например:

• Метод бисекции — отличается гарантированной сходимостью, однако с довольно низкой скоростью
• Метод хорд — не требует вычисления производной функции, обладает неплохой скоростью сходимости

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

4. Определяем параметры траектории пушечного ядра

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

Делаем это исходя из физического смысла задачи — после сразу выстрела высота полета снаряда неотрицательна. Перебираем все моменты времени, начиная от нуля, до тех пор, пока высота не станет отрицательна. Перебираем с достаточно крупным шагом (1 секунда) чтобы процедура не была слишком длительной, ведь на каждом шаге мы заново интегрируем дифференциальные уравнения движения, что весьма накладно сточки зрения вычислительных затрат. Как только высота станет отрицательной, заканчиваем перебор. Корень уравнения h(t) = 0 находится где-то внутри интервала [a, b]. Начальное приближение берем в середине этого интервала

Теперь отдаем уравнение на съедение процедуре решения нелинейных уравнений Octave

Функция fsolve() на вход принимает функцию, описывающую левую часть уравнения f(x) = 0 и значение начального приближения. Возвращает значение корня, вычисленное с заданной точностью. С какой точностью? Пока не будем задаваться этим вопросом и воспользуемся настройками по-умолчанию, на данном этапе они нас устраивают.

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

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

В командном окне можно увидеть результаты работы программы

а также посмотреть, с какой точностью были решены уравнения

Для наших учебных целей точность вполне приемлема.

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

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

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

Дано нелинейное уравнение:

( 4.1)

Необходимо решить это уравнение, т. е. найти его корень .

Если функция имеет вид многочлена степени m,

где a_i — коэффициенты многочлена, , то уравнение f(x)=0 имеет m корней (рис. 4.2).

Если функция f(x) включает в себя тригонометрические или экспоненциальные функции от некоторого аргумента x , то уравнение (4.1) называется трансцендентным уравнением .

Такие уравнения обычно имеют бесконечное множество решений.

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

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

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

Большинство употребляющихся приближенных методов решения уравнений являются, по существу, способами уточнения корней. Для их применения необходимо знание интервала изоляции [a,b] , в котором лежит уточняемый корень уравнения (рис. 4.3).

Процесс определения интервала изоляции [a,b] , содержащего только один из корней уравнения, называется отделением этого корня.

Процесс отделения корней проводят исходя из физического смысла прикладной задачи, графически, с помощью таблиц значений функции f(x) или при помощи специальной программы отделения корней. Процедура отделения корней основана на известном свойстве непрерывных функций: если функция непрерывна на замкнутом интервале [a,b] и на его концах имеет различные знаки, т.е. f(a)f(b) , то между точками a и b имеется хотя бы один корень уравнения (1). Если при этом знак функции f'(x) на отрезке [a,b] не меняется, то корень является единственным на этом отрезке.

Процесс определения корней алгебраических и трансцендентных уравнений состоит из 2 этапов:

1. отделение корней, — т.е. определение интервалов изоляции [a,b] , внутри которого лежит каждый корень уравнения;
2. уточнение корней, — т.е. сужение интервала [a,b] до величины равной заданной степени точности .

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

Решение нелинейных уравнений

Уравнения, в которых содержатся неизвестные функции, произведенные в степень больше единицы, называются нелинейными.
Например, y=ax+b – линейное уравнение, х^3 – 0,2x^2 + 0,5x + 1,5 = 0 – нелинейное (в общем виде записывается как F(x)=0).

Системой нелинейных уравнений считается одновременное решение нескольких нелинейных уравнений с одной или несколькими переменными.

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

Рассмотрим несколько методов уточнения корней с определенно заданной точностью.

Методы численного решения нелинейных уравнений

Метод половинного деления.

Суть метода половинного деления заключается в делении интервала [a,b] пополам (с=(a+b)/2) и отбрасывании той части интервала, в которой отсутствует корень, т.е. условие F(a)xF(b)

Рис.1. Использование метода половинного деления при решении нелинейных уравнений.

Рассмотрим пример. Необходимо решить уравнение х^3 – 0,2x^2 + 0,5x + 1,5 = 0 с точностью до e 0, то начала отрезка a переносится в x (a=x), иначе, конец отрезка b переносится в точку x (b=x). Полученный отрезок делим опять пополам и т.д. Весь произведенный расчет отражен ниже в таблице.

Рис.2. Таблица результатов вычислений

В результате вычислений получаем значение с учетом требуемой точности, равной x=-0,946

При использовании метода хорд, задается отрезок [a,b], в котором есть только один корень с установленной точностью e. Через точки в отрезке a и b, которые имеют координаты (x(F(a);y(F(b)), проводится линия (хорда). Далее определяются точки пересечения этой линии с осью абсцисс (точка z).
Если F(a)xF(z)

Рис.3. Использование метода хорд при решении нелинейных уравнений.

Рассмотрим пример. Необходимо решить уравнение х^3 – 0,2x^2 + 0,5x + 1,5 = 0 с точностью до e 0;

Определим вторую производную F’’(x) = 6x-0,4.

F’’(-1)=-6,4 0 соблюдается, поэтому для определения корня уравнения воспользуемся формулой:

, где x0=b, F(a)=F(-1)=-0,2

Весь произведенный расчет отражен ниже в таблице.

Рис.4. Таблица результатов вычислений

В результате вычислений получаем значение с учетом требуемой точности, равной x=-0,946

Метод касательных (Ньютона)

Данный метод основывается на построении касательных к графику, которые проводятся на одном из концов интервала [a,b]. В точке пересечения с осью X (z1) строится новая касательная. Данная процедура продолжается до тех пор, пока полученное значение не будет сравним с нужным параметром точности e (F(zi)

Рис.5. Использование метода касательных (Ньютона) при решении нелинейных уравнений.

Рассмотрим пример. Необходимо решить уравнение х^3 – 0,2x^2 + 0,5x + 1,5 = 0 с точностью до e 0 выполняется, поэтому расчеты производим по формуле:

Весь произведенный расчет отражен ниже в таблице.

Рис.6. Таблица результатов вычислений

В результате вычислений получаем значение с учетом требуемой точности, равной x=-0,946

Если материал был полезен, вы можете отправить донат или поделиться данным материалом в социальных сетях: