Решение систем дифференциальных уравнений программирование

Как написать программу решения дифференциального уравнения… ( C++ )

Для численного решения обыкновенных дифференциальных уравнений различают задачи с начальными условиями (ЗНУ) и граничными условиями (ЗГУ).

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

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

Для численного решения обыкновенных дифференциальных уравнений применяют метод Рунге-Кутта, с наиболее часто используемым 4-ым порядком точности. Давайте рассмотрим программную реализацию данного метода Рунге-Кутта, и Вы увидите, что написать программу для решения дифференциального уравнения не составляет особого труда.

Вариант решения задачи рассмотрю на следующем примере:

Условия задачи:
Пусть выстрел из орудия произведен с начальной скоростью V₀, под углом к горизонту α, с высоты Н₀ расположения орудия, т.к. в реальности огонь может вестись с холма или из капонира (т.е. ниже уровня земли).
Считаем, что снаряд имеет форму шара с радиусом r, изготовлен из материала, имеющего определенную плотность ρ.
Построить траекторию полета снаряда Y(x) ,
указать максимальную высоту полета H_k , дальность падения снаряда X_k и время полета t_k , построить график скорости V(t) на отрезке [0,t_k].

Таким образом, исходные данные, которые пользователь может задать на форме:
Начальная скорость V₀, м/с2
Начальная высота H₀, м
Угол выстрела α, ° (град)
Плотность материала ρ, кг/м3
Радиус r, м

При построении математической модели условимся, что ось Оx системы координат направлена горизонтально в направлении выстрела, а ось Oy — вертикально вверх. Вектор скорости снаряда V(t) за время полета будет изменяться как по величине, так и по направлению, поэтому в модели рассматриваем его проекции на координатные оси. Горизонтальную составляющую скорости в момент времени t обозначим V_x(t), а вертикальную – V_y(t).

Пусть поверхность Земли плоская. Согласно законам механики, при сделанных предположениях движения тела в горизонтальном направлении является равномерным, а в вертикальном – равнозамедленным или равноускоренным с ускорением свободного падения g.

Если с силой тяжести F_T все достаточно просто (она свой вектор не меняет ни по величине, ни по направлению), то сила лобового сопротивления F_C , действующая на снаряд, пропорциональна квадрату скорости движения тела. Обозначим через F_X и F_Y горизонтальную и вертикальную проекции вектора силы лобового сопротивления,
причем F_X/F= V_X/V, F_Y/F= V_Y/V.

Значение силы лобового сопротивления F= -b·V² (пропорционально квадрату скорости тела). Коэффициент b=0.5·C·S·ρ, где C – коэффициент лобового сопротивления (для многих задач баллистики C≈0.15), S – площадь поперечного сечения (S=πr²), ρ — плотность воздуха (ρ=1,29 кг/м3).

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

Метод Рунге-Кутта предполагает многократное вычисление значения производной искомой функции по имеющейся формуле (из уравнения), поэтому имеет смысл …

Код функций будет выглядеть так:

// функция Нахождение горизонтальной проекции скорости
//по первому уравнению системы
double Form1::fvx( double vy, double vx )
<
return -b*vx*sqrt(vx*vx+vy*vy) / m;
>

// функция Нахождение вертикальной проекции скорости
//по второму уравнению системы
double Form1::fvy( double vy, double vx )
<
return -b*vy*sqrt(vx*vx+vy*vy) / m — g;
>

Шаг интегрирования у меня задается на форме.
Сейчас нам предстоит вычислить значения нескольких функций (V_x, V_y, V ) в точках интервала с шагом. В моем примере у интервала задано начало х=0, а конечная точка интервала будет определена в процессе вычисления ( высота полета ядра стала Void Form1::Runge_Kutta(void)
<
double k1,k2,k3,k4, l1,l2,l3,l4;

pY[0]=H; pX[0]=0; pt[0]=0; //массивы-координаты: высота, дальность и время
pVx[0]=Vx; pVy[0]=Vy; pV[0]=V; //массивы- скорости: проекции на
//горизонталь и вертикаль и полная скорость (величина вектора)
bool vzbool=true;//взлет
int i=1;

//расчет по модели и заполнение массивов
while( (pY[i-1]>-0.00001 || vzbool) && i pY[i-1]) iMax=i; //сохранение номера узла с максимальной высотой
else vzbool=false;//падение

i++;
>
n=i-1; //количество реальных шагов
>

где:
int iMax; //узел с макс.высотой полета

double b; //коэф.пропорциональности Силы лобового сопротивления
double m; //масса ядра
double H; //уровень расположения орудия в момент выстрела
double V,Vx,Vy; //начальная скорость и ее проекции на оси

В результате работы этой подпрограммы произойдет численное решение задачи Коши для системы обыкновенных дифференциальных уравнений и будут получены значения функций V_x(t_i), V_y(t_i) в точках t_i=i·h, i=1,2. ; h – шаг метода.

Как видим, после получения нового значения скорости V_x(t_i)
рассчитывается координата X(t_i)=X(t_i-1)+h·V_x(t_i), где h= t_i-t_i-1 = const.
Кроме того, параллельно рассчитывается значение высоты Y(t_i)=Y(t_i-1)+h·V_y(t_i),
где h= t_i-t_i-1 = const по найденным значениям V_y(t_i).
Когда будет получено значение Y(t_i)

Символьное решение линейных дифференциальных уравнений и систем методом преобразований Лапласа c применением SymPy

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

В данной публикации предлагаю рассмотреть функции прямого и обратного преобразования Лапласа из библиотеки SymPy, которые позволяют использовать метод Лапласа для решения дифференциальных уравнений и систем средствами Python.

Сам метод Лапласа и его преимущества при решении линейных дифференциальных уравнений и систем широко освещены в литературе, например в популярном издании [1]. В книге метод Лапласа приведен для реализации в лицензионных программных пакетах Mathematica, Maple и MATLAB (что подразумевает приобретение учебным заведением этого ПО) на выбранных автором отдельных примерах.

Попробуем сегодня рассмотреть не отдельный пример решения учебной задачи средствами Python, а общий метод решения линейных дифференциальных уравнений и систем с использованием функций прямого и обратного преобразования Лапласа. При этом сохраним обучающий момент: левая часть линейного дифференциального уравнения с условиями Коши будет формироваться самим студентом, а рутинная часть задачи, состоящая в прямом преобразовании Лапласа правой части уравнения, будет выполняться при помощи функции laplace_transform().

История об авторстве преобразований Лапласа

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

В данном случае следующим после Эйлера был французский математик Пьер Симон де Лаплас (Pierre Simon de Laplace (1749-1827)). Именно он использовал такие интегралы в своей работе по теории вероятностей. Самим Лапласом не применялись так называемые «операционные методы» для нахождения решений дифференциальных уравнений, основанные на преобразованиях Лапласа (изображениях по Лапласу). Эти методы в действительности были обнаружены и популяризировались инженерами-практиками, особенно английским инженером-электриком Оливером Хевисайдом (1850-1925). Задолго до того, как была строго доказана справедливость этих методов, операционное исчисление успешно и широко применялось, хотя его законность ставилось в значительной мере под сомнение даже в начале XX столетия, и по этой теме велись весьма ожесточенные дебаты.

Функции прямого и обратного преобразования Лапласа

Эта функция возвращает (F, a, cond), где F(s) есть преобразование Лапласа функции f(t), a Текст программы

Время на обратное визуальное преобразование Лапласа: 2.68 s

Обратное преобразование Лапласа часто используется при синтезе САУ, где Python может заменить дорогостоящих программных “монстров” типа MathCAD, поэтому приведенное использование обратного преобразования имеет практическое значение.

Преобразование Лапласа от производных высших порядков для решения задачи Коши

Если a и b — константы, то

для всех s, таких, что существуют оба преобразования Лапласа (изображения по Лапласу) функций f(t) и q(t).

Проверим линейность прямого и обратного преобразований Лапласа с помощью ранее рассмотренных функций laplace_transform() и inverse_laplace_transform(). Для этого в качестве примера примем f(t)=sin(3t), q(t)=cos(7t), a=5, b=7 и используем следующую программу.

(7*s**3 + 15*s**2 + 63*s + 735)/((s**2 + 9)*(s**2 + 49))
(7*s**3 + 15*s**2 + 63*s + 735)/((s**2 + 9)*(s**2 + 49))
True
5*sin(3*t) + 7*cos(7*t)
5*sin(3*t) + 7*cos(7*t)

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

Если предположить, что удовлетворяет условиям первой теоремы, то из этой теоремы будет следовать, что:

Повторение этого вычисления дает

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

Применяя соотношение (3), содержащее преобразованные по Лапласу производные искомой функции с начальными условиями, к уравнению (1), можно получить его решение по методу, специально разработанному на нашей кафедре при активной поддержке Scorobey для библиотеки SymPy.

Метод решения линейных дифференциальных уравнений и систем уравнений, основанный на преобразованиях Лапласа, с использованием библиотеки SymPy

где — приведенное начальное положение массы, — приведенная начальная скорость массы.

Упрощённая физическая модель, заданная уравнением (4) при ненулевых начальных условиях [1]:

Система, состоящая из материальной точки заданной массы, закрепленной на пружине, удовлетворяет задаче Коши (задаче с начальными условиями). Материальная точка заданной массы первоначально находится в покое в положении ее равновесия.

Для решения этого и других линейных дифференциальных уравнений методом преобразований Лапласа удобно пользоваться следующей системой, полученной из соотношений (3):

Последовательность решения средствами SymPy следующая:

загружаем необходимые модули и явно определяем символьные переменные:

указываем версию библиотеки sympy, чтобы учесть ее особенности. Для этого нужно ввести такие строки:

по физическому смыслу задачи переменная времени определяется для области, включающей ноль и положительные числа. Задаём начальные условия и функцию в правой части уравнения (4) с её последующим преобразование по Лапласу. Для начальных условий необходимо использовать функцию Rational, поскольку использование десятичного округления приводит к ошибке.

пользуясь (5), переписываем преобразованные по Лапласу производные, входящие в левую часть уравнения (4), формируя из них левую часть этого уравнения, и сравниваем результат с правой его частью:

решаем полученное алгебраическое уравнение относительно преобразования X(s) и выполняем обратное преобразование Лапласа:

осуществляем переход из работы в библиотеке SymPyв библиотеку NumPy:

строим график обычным для Python методом:

Получаем:
Версия библиотеки sympy – 1.3

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

При этом учебное значение метода решения сохраняется за счёт необходимости использования системы (5) и перехода в NumPy для исследования решения более производительными методами.

Для дальнейшей демонстрации метода решим систему дифференциальных уравнений:

с начальными условиями

Упрощённая физическая модель, заданная системой уравнений (6) при нулевых начальных условиях:

Таким образом, сила f(t) внезапно прилагается ко второй материальной точке заданной массы в момент времени t = 0, когда система находится в покое в ее положении равновесия.

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

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

Рассмотрим решение линейного дифференциального уравнения четвёртого порядка с нулевыми начальными условиями:

Решим линейное дифференциальное уравнение четвёртого порядка:

с начальными условиями , , .

Функции для решения ОДУ

Для имеющих аналитическое решение ОДУ и систем ОДУ применяется функция dsolve():
sympy.solvers.ode.dsolve(eq, func=None, hint=’default’, simplify=True, ics=None, xi=None, eta=None, x0=0, n=6, **kwargs)

Давайте сравним производительность функции dsolve() с методом Лапласа. Для примера возьмём следующее дифференциальное уравнение четвёртой степени с нулевыми начальными условиями:

Время решения уравнения с использованием функции dsolve(): 1.437 s

Время решения уравнения с использованием преобразования Лапласа: 3.274 s

Итак, функция dsolve() (1.437 s) решает уравнение четвёртого порядка быстрее, чем выполняется решение по методу преобразований Лапласа (3.274 s) более чем в два раза. Однако при этом следует отметить, что функция dsolve() не решает системы дифференциальных уравнений второго порядка, например, при решении системы (6) с использованием функция dsolve() возникает ошибка:

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

Для того чтобы найти необходимый метод решения дифференциальных уравнений с помощью функции dsolve(), нужно использовать classify_ode(eq, f(x)), например:

Eq(f(x), C1*sin(x) + C2*cos(x))
(‘nth_linear_constant_coeff_homogeneous’, ‘2nd_power_series_ordinary’)
(‘separable’, ‘1st_exact’, ‘almost_linear’, ‘1st_power_series’, ‘lie_group’, ‘separable_Integral’, ‘1st_exact_Integral’, ‘almost_linear_Integral’)
[Eq(f(x), -acos((C1 + Integral(0, x))*exp(-Integral(-tan(x), x))) + 2*pi), Eq(f(x), acos((C1 + Integral(0,x))*exp(-Integral(-tan(x), x))))]

Таким образом, для уравнения eq=Eq(f(x).diff(x,x)+f(x),0) работает любой метод из первого списка:

Для уравнения eq = sin(x)*cos(f(x)) + cos(x)*sin(f(x))*f(x).diff(x) работает любой метод из второго списка:

separable, 1st_exact, almost_linear,
1st_power_series, lie_group, separable_Integral,
1st_exact_Integral, almost_linear_Integral

Чтобы использовать выбранный метод, запись функции dsolve() примет вид, к примеру:

Вывод:

Данная статья ставила своей целью показать, как использовать средства библиотек SciPy и NumPy на примере решения систем линейных ОДУ операторным методом. Таким образом, были рассмотрены методы символьного решения линейных дифференциальных уравнений и систем уравнений методом Лапласа. Проведен анализ производительности этого метода и методов, реализованных в функции dsolve().

1. Дифференциальные уравнения и краевые задачи: моделирование и вычисление с помощью Mathematica, Maple и MATLAB. 3-е издание.: Пер. с англ. — М.: ООО «И.Д. Вильяме», 2008. — 1104 с.: ил. — Парал. тит. англ.
2. Использование обратного преобразования Лапласа для анализа динамических звеньев систем управления

Решение систем обыкновенных дифференциальных уравнений в среде MATLAB. Часть 1

В среде MATLAB можно решать системы диффуров с начальными условиями, краевые задачи, а также решать дифференциальные уравнения в частных производных с помощью инструмента PDE toolbox.

В данном обзоре речь пойдет лишь о системах дифференциальных уравнений с начальными условиями, то есть о задаче Коши. В англоязычной литературе это называется Initial Value Problem.

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

Решать системы обыкновенных дифференциальных уравнений можно как в MATLAB, так и в Simulink.

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

Выбор ваш должен зависеть от задачи. Если Вы, например, хотите смоделировать какой-либо объект управления, описываемый системой диффуров, то в данном случае имеет смысл использовать именно Simulink, так как Вам, впоследствии, понадобиться синтез, например, системы управления, и Simulink подойдет здесь как нельзя лучше.

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

Рассмотрим синтаксис решателей matlab.В качестве аргументов следует подать правую часть системы в виде MATLAB-функции.

На рисунке показан требуемый вид системы, когда выражены старшие производные.

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

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

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

Далее задаются начальные условия. Значения всех неизвестных искомых переменных в начале расчёта задаются в виде столбца соответствующей размерности.

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

Вы могли заметить, что название функции — odeXY – это обозначение для всех решателей, которых всего 8 штук. В данном ролике мы пользоваться решателем ode45, соответствующего численному по методу Дормана-Принса 4(5). Этого решателя достаточно для подавляющего большинства задач. Остальные решатели будут подробно рассмотрены в приложении к задачам соответствующих типов позже.

Перейдем к примерам.

Рассмотрим 2 примера:

• решение дифференциального уравнения первого порядка.
• решение системы двух дифференциальных уравнений второго порядка.

В качестве уравнение первого порядка рассмотрим логистическое уравнение Ферхюльста, которое описывает динамику численности популяции. Суть уравнения такова: скорость прироста населения y пропорциональна количеству населения, однако лимитирована максимальной численностью популяции.

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

Решение этого дифференциального уравнения выглядит следующим образом.

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

Перейдем в окно MATLABа и посмотрим, как это выглядит.

Так выглядит скрипт:

Так выглядит график решения дифференциального уравнения:

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

Уравнения показаны на рисунке. Но вид системы отличается от требуемого, в том числе потому, что в нём присутствуют вторые производные. Для приведения системы в требуемый вид выполним 2 простых шага:

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

Настало время её записать.

Итак, мы имеем систему, параметры, интервал времени и начальные условия. Решим же эту задачу скорее.

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

Эту функцию можно располагать как в самом скрипте решения в самом его конце, так и в виде отдельного m-файла.

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

Теперь рассмотрим скрипт самого решения.

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

Сообразно с уже разобранным ранее синтаксисом укажем функцию pendulum_np, интервал времени и начальные условия.

Перейдем теперь в окно MATLAB и посмотрим решение.

Так выглядит скрипт:

Запускаем скрипт и получаем графики:

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

На рисунке показана функция MATLAB, которая соответствует движению подвешенного на пружине шара, однако можно заметить, что эта функция теперь имеет на 5 аргументов больше.

Параметры задаются в скрипте, а при вызове функции мы обращаемся к уже известному оператору-собаке, которая превращает функцию семи переменных pendulum_n в функцию двух переменных t и X. Вот и всё.

Я вам очень рекомендую разобраться с тем, как работает оператор-собака. В хелпе он называется function-handle. Разобравшись с ним Вам будет работать в среде MATLAB ещё проще и ещё приятнее.

Вывод: не так страшно решать диффуры

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

Их можно, с одной стороны, разделить по степени жёсткости, а с другой стороны, по структуре самой системы.

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

Видеообзор по теме решения систем Д/У доступен по ссылке.