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

Введение

Содержание цикла
Введение
  1. Численные методы решения ОДУ
  2. Введение в GNU Octave
  3. Задача внешней баллистики


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

mjscamdvtcf7uapy793tego4wy4.jpeg

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

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

$\begin{align} &x = x(t) \\ &y = y(t) \\ &v_x = v_x(t) \\ &v_y = v_y(t) \end{align}$


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

$y(t) = 0$


относительно времени, то мы найдем момент времени $t^{*}$ в который снаряд упал на землю. Координата $x(t^{*})$ и будет дальность полета снаряда.

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

ur_te3mdap_g5-gq4hzct1hgfb4.png

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

$\dot y(t) = 0$


или

$v_y(t) = 0$


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

$h_{\max} = y(t^{**})$

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

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

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

$x + 2 = e^{x} $


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

На самом деле ответ есть
Проведем некоторые эквивалентные преобразования

$ \begin{align} &x + 2 = e^{x} \\ &(x + 2) \, e^{-x} = 1 \\ \end{align} $


Введем замену $u = x + 2$, тогда
\begin{align}
&u \, e^{2 — u} = 1 \\
&u \, e^{-u} = e^{-2} \\
&-u \, e^{-u} = -e^{-2}
\end{align}
Пусть теперь $w = -u$, тогда

$w \, e^{w} = -e^{-2}$


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

$z = w \, e^{w}$


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

$w = W(z)$


Не в даваясь в теорию ФКП, в которой я мало смыслю, скажу, что число $-e^{-2}$ попадает в интервал $[-e^{-1}; 0)$ в котором функция Ламберта многозначна, значит корня будет два

$w_1 = W_0(-e^{-2}), \quad w_2 = W_{-1}(-e^{-2})$


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

$x_1 = -W_0(-e^{-2}) - 2, \quad x_2 = -W_{-1}(-e^{-2}) - 2$


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

$x_1 = -1.841405660, \quad x_2 = 1.146193221$



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

Графическое решение уравнения
nxkkwlsrpncp4hx575zmywmokam.png


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

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

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

%-------------------------------------------------------------------------------
% Параметры снаряда
%-------------------------------------------------------------------------------

% плотность материала снаряда
global gamma = 7800; 
% коэффициент сопротивления формы (для шара)
global c_f = 0.47;
% калибр
global d = 0.1;

%-------------------------------------------------------------------------------
% Параметры стрельбы
%-------------------------------------------------------------------------------

% Начальная скорость
global v0 = 400;
% Угол наклона ствола пушки (в градусах!)
global alpha = 35;


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

%-------------------------------------------------------------------------------
%   Расчет точки фазовой траектории для произвольного момента времени
%-------------------------------------------------------------------------------
function Y = MotionOdeSolve(v0, alpha, t)
   
  % Начальные условия
  x0 = 0;
  y0 = 0;
  vx0 = v0 * cos(deg2rad(alpha));
  vy0 = v0 * sin(deg2rad(alpha));
  
  Y0 = [x0; y0; vx0; vy0];
  
  % Решаем ОДУ движения
  solv = lsode("f_air", Y0, [0; t]);
  
  % Формируем вектор фазовых координат для момента времени t
  Y = [solv(2,1); solv(2,2); solv(2,3); solv(2,4)];
  
endfunction


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

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

%-------------------------------------------------------------------------------
% Вычисление высоты полета снаряда для произвольного момента времени
%-------------------------------------------------------------------------------
function h = height(t)
  
  global v0;
  global alpha;
  
  Y = MotionOdeSolve(v0, alpha, t);
  
  h = Y(2);
  
endfunction


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

t1 = 10;
t2 = 20;
printf("h(%f) = %f, h(%f) = %f\n", t1, height(t1), t2, height(t2));


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

>> cannon
h(10.000000) = 855.013452, h(20.000000) = 500.106649


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

%-------------------------------------------------------------------------------
% Горизонтальное удаление снаряда от позиции стрельбы для произвольного t
%-------------------------------------------------------------------------------
function d = dist(t)
  
  global v0;
  global alpha;
  
  Y = MotionOdeSolve(v0, alpha, t);
  
  d = Y(1);
  
endfunction


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

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

$f(x) = 0$


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

$x + 2 = e^{x}$


превращается в 

$x + 2 - e^{x} = 0$


где $f(x) = x + 2 - e^{x}$. Корни этого, эквивалентного уравнения, равны корням исходного. Если мы построим график функции f (x), то увидим такую картинку

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

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

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


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

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

$x = \varphi(x)$


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

$x_{i+1} = \varphi(x_i)$


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

$x_{i+1} = e^{x_i} - 2$


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

$\begin{align} &x_1 = e^{x_0} - 2 = e^{-2} - 2 = -1.8647 \\ &x_2 = e^{x_1} - 2 = e^{-1.8647} - 2 = -1.8451 \\ &x_3 = e^{x_2} - 2 = e^{-1.8451} - 2 = -1.8412 \\ &x_4 = e^{x_3} - 2 = e^{-1.8412} - 2 = -1.8415 \\ &x_5 = e^{x_4} - 2 = e^{-1.8415} - 2 = -1.8414 \\ &x_6 = e^{x_5} - 2 = e^{-1.8414} - 2 = -1.8414 \end{align}$


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

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

$f(x) = f(x_0) + f'(x_0) \, (x - x_0) + \frac{1}{2} \, f''(x_0) \, (x - x_0)^2 + ...$


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

$f(x) \approx f(x_0) + f'(x_0) \, (x - x_0) $


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

$x_1 = x_0 - \frac{f(x_0)}{f'(x_0)}$


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

$x_{i+1} = x_i - \frac{f(x_i)}{f'(x_i)}$


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

$\begin{align} &f(x) = x + 2 - e^{x} \\ &f'(x) = 1 - e^{x} \end{align}$


итерационная формула

$x_{i+1} = x_{i} - \frac{x_i + 2 - e^{x_i}}{1 - e^{x_i}}$


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

$\begin{align} &x_1 = x_0 - \frac{x_0 + 2 - e^{x_0}}{1 - e^{x_0}} = 1 - \frac{1 + 2 - e}{1 - e} = 1.1640 \\ &x_2 = x_1 - \frac{x_1 + 2 - e^{x_1}}{1 - e^{x_1}} = 1.1640 - \frac{1.1640 + 2 - e^{1.1640}}{1 - e^{1.1640}} = 1.1464 \\ &x_3 = x_2 - \frac{x_2 + 2 - e^{x_2}}{1 - e^{x_2}} = 1.1464 - \frac{1.1464 + 2 - e^{1.1464}}{1 - e^{1.1464}} = 1.1462 \\ &x_4 = x_3 - \frac{x_3 + 2 - e^{x_3}}{1 - e^{x_3}} = 1.1462 - \frac{1.1462 + 2 - e^{1.1462}}{1 - e^{1.1462}} = 1.1462 \end{align}$


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

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

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


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

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

%
% Решаем задачу поиска дальности полета снаряда. Следует решить уравнение
%   height(t) = 0
% относительно времени. 
%

% Поиск интервала изоляции корня
t0 = 0;
deltaT = 1.0;

t = t0;

while (height(t) >= 0)
  t = t + deltaT;  
endwhile

b = t;
a = t - deltaT;

% Первое приближение корня
t1 = t - deltaT / 2;


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

$t_1 = \frac{a + b}{2} = t - \frac{\Delta t}{2}$


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

% Уточняем корень
t_end = fsolve("height", t1);


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

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

% Определяем дистанцию стрельбы
d_max = dist(t_end);
% выводим её в терминал
printf("Shot distance: %f, m\n", d_max);


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

%-------------------------------------------------------------------------------
% Вычисление вертикальной проекции вектора скорости
%-------------------------------------------------------------------------------
function v = vy(t)
  
  global v0;
  global alpha;
  
  Y = MotionOdeSolve(v0, alpha, t);
  
  v = Y(4);
  
endfunction

%
% Решаем задачу поиска максимальной высоты подъема снаряда. Следует решить 
% уравнение
%
% vy(t) = 0
%

% Ищем интервал изоляции корня
t = t0;

while (vy(t) >= 0)
  t = t + deltaT;  
endwhile

% Первое приближение корня
t1 = t - deltaT / 2;

% Уточняем корень
t_hmax = fsolve("vy", t1);

% Находим максимальную высоту подъема снаряда
h_max = height(t_hmax);

printf("Maximal height: %f, m\n", h_max);


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

>> cannon
h(10.000000) = 855.013452, h(20.000000) = 500.106649
Shot distance: 2840.991508, m
Maximal height: 857.963929, m


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

>> height(t_end)
ans =   -1.2328e-06
>> vy(t_hmax)
ans =    4.5117e-09


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

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

© Habrahabr.ru