Решение дифференциальных уравнений

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

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

Общее решение обыкновенного дифференциального уравнения порядка имеет вид:

где , , ... , --- произвольные постоянные (постоянные интегрирования). Каждый частный выбор этих постоянных дает частное решение. В задаче Коши (начальной задаче) требуется найти частное решение, удовлетворяющее начальным условиям:

по которым определяются постоянных , , ... ,

Вспомним как решаются дифференциальные уравнения на примере при начальных условиях: при .

Определим тип уравнения : дифференциальное уравнение с разделяющимися переменными. Поэтому для решения переносим переменные:

далее интегрируем:

общее решение выглядит следующим образом:

чтобы определить константу интегрирования поставим начальные условия в уравнение:

Частное решение:

избавимся от логарифма и выразим функцию в явном виде:

Численное решение дифференциальных уравнений

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

using DifferentialEquations;

[ Info: Precompiling PreallocationToolsSparseConnectivityTracerExt [f3463b80-2822-5513-b06b-5bb94f50d0b8] 
[ Info: Precompiling ArrayInterfaceBandedMatricesExt [26e938bc-0cd5-5679-9003-44616cbf91d3]

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

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

f(y, p, x) = -1 / 3 * y

f (generic function with 1 method)

Для решения также необходимо задать начальное условие

y₀ = 10.0

10.0

и диапазон на котором будем решать дифференциальное уравнение

x = (0.0, 10.0)

(0.0, 10.0)

Далее формируется проблема в виде задачи Коши:

prob = ODEProblem(f, y₀, x)

ODEProblem with uType Float64 and tType Float64. In-place: false
Non-trivial mass matrix: false
timespan: (0.0, 10.0)
u0: 10.0

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

sol = solve(prob);

Количество значений функции в массиве определяется методом решения уравнения и заданной точностью.
Значения функций хранятся в sol.u

sol.u

10-element Vector{Float64}:
 10.0
  9.593251647361855
  7.961188151479715
  5.9540975559518445
  4.27692687235221
  2.835132447794858
  1.7880393651249318
  1.0627000769015806
  0.604618513074874
  0.35674754779955575

соответствующие значений хранятся в sol.t

sol.t

10-element Vector{Float64}:
  0.0
  0.12457558532377848
  0.6840205291300432
  1.5555165873605028
  2.5480519184065056
  3.781492790133791
  5.164405599413932
  6.72533998547542
  8.417276529952703
 10.0

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

sol(5)

1.8887602150416114

using Plots

y(x) = 10 * exp(-1/3 * x)   #аналитическое решение

y (generic function with 1 method)

plot(sol.t, sol.u, label="численное решение")
plot!(sol.t, y.(sol.t), label = "аналитическое решение")

Решение систем уравнений

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

Пример, решим систему уравнений

Решить систему дифференциальных уравнений в интервале [0; 10] с начальными условиями и

Результат решения представить графически.

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

Функции записываются в виде массива, y - первая функция (u[1]), z - вторая функция (u[2]), x - переменная по которой ведется интегрирование (t).

function fsys(du, u, p, t)
	du[1] = sin(u[1] - u[2]) * t
	du[2] = t  / u[1] #sin(u[1] - u[2]) * t 
end

fsys (generic function with 1 method)

tspan = (0,10)

(0, 10)

u0 = [1, -1]

2-element Vector{Int64}:
  1
 -1

sysprob = ODEProblem(fsys, u0, tspan);

syssol = solve(sysprob);

plot(syssol)

Дифференциальные уравнения высших порядков

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

На диапазоне от 0 до 10, при начальных условиях и

Сделаем замену , и

function fsys₂(du, u, p, t)
	du[1] = 0.1 * u[2] - u[1] - t
	du[2] = u[2]
end

fsys₂ (generic function with 1 method)

tspan₂ = (0,10)

(0, 10)

u0₂ = [1, -1]

2-element Vector{Int64}:
  1
 -1

sysprob₂ = ODEProblem(fsys₂, u0₂, tspan₂);

syssol₂ = solve(sysprob₂);

y₂(x) = syssol₂(x)[1]

y₂ (generic function with 1 method)

x₂ = 0:0.2:10

0.0:0.2:10.0

plot(x₂, y₂.(x₂), label="y(x)")

Краевая задача

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

Пример, решим систему дифференциальных уравнений

на диапазоне от 0 до 10 при граничных условиях и

using BoundaryValueDiffEq

Ситема уравнений задается аналогичным образом

function fsys₃(du, u, p, t)
	du[1] = u[1] - u[2]
	du[2] = t - u[1]
end

fsys₃ (generic function with 1 method)

tspan₃ = (0,10)

(0, 10)

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

преобразуется к виду:

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

function bc₃(r, u, du, t)
	r[1] = u(0)[1] - 1
	r[2] = u(10)[2] - 5
end

bc₃ (generic function with 1 method)

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

uᵢₙ = [1.0, 1.0]

2-element Vector{Float64}:
 1.0
 1.0

bvp1 = BoundaryValueDiffEq.BVProblem(fsys₃, bc₃, uᵢₙ, tspan₃);

sol₃ = solve(bvp1, BoundaryValueDiffEq.MIRK4(); dt = 0.05);

plot(sol₃)

Пример решения задачи из механики

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

По закону Ньютона:

Рассмотрим пример решения задачи маятника.

В декартовой системе координат направим ось y вниз. Система уравнений будет иметь вид

где N - сила натяжения нити, L - длина нити.
Эти уравнения необходимо дополнить уравнением описывающим, что нить имеет постоянную длину: .

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

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

Запишем данную систему в виде функции (значения g и L возьмем в качестве параметров):

function mf(du, u, p, t)
	g = p[1]
	L = p[2]
	#
	x = u[1]
	υₓ = u[2]
	y = u[3]
	υᵧ = u[4]   #нижнего индекса y нет в юникоде, поэтому использую похужую в написании гамму

	du[1] = υₓ
	du[2] = - g  / L^2 * x * y - 1 / L^2 * x * (υₓ^2 + υᵧ^2)
	du[3] = υᵧ
	du[4] = (- g / L^2 * y^2 - 1 / L^2 * y * (υₓ^2 + υᵧ^2) + g)
end

mf (generic function with 1 method)

Задаем длину нити:

L = 1.0

1.0

Задаем начальные угол на который отклонится груз

θ = -10 / 180 * pi

-0.17453292519943295

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

init = [cos(θ) * L, 0.0, sin(θ) * L, 0.0]

4-element Vector{Float64}:
  0.984807753012208
  0.0
 -0.17364817766693033
  0.0

Решим уравнение в диапазоне до 3 сек

τ = (0,3)

(0, 3)

mprob = ODEProblem(mf, init, τ, [9.8, L]);

msol = solve(mprob);

Из решения возьмем необходимые значения при времени с шагом в 0.01 сек.

Тут также мы перевернем ось y т.к. при записи уравнений мы ее направили вниз.

τᵣ = 0:0.01:3;

xᵣ = collect(msol(τ)[1] for τ in τᵣ);  #первое значение из решения соотвествует x

υₓ = collect(msol(τ)[2] for τ in τᵣ);  #второе скорости по x и т.д.

yᵣ = -collect(msol(τ)[3] for τ in τᵣ);

υᵧ = -collect(msol(τ)[4] for τ in τᵣ);

Нарисуем траекторию движения шара.

scatter(xᵣ, yᵣ, xlims= (-1.05,1.05), ylims = (-1.05,1.05), label = "", aspect_ratio =1)
# xlims и ylims задают диапазон, в котором будет строится график (по умолчанию масштаб выберется автоматически)
# aspect_ratio =1 означает, что масштабы по x и y будут одинаковыми. Это необходимо чтобы адекватно воспринималась длина нити.

Но данная картинка не очень интересна. Лучше построим изменение траектории по времени. Для этого воспользуемся возможностью создавать анимации в пакете Plots.

anim = @animate for i in 1:length(τᵣ)
	scatter([0,xᵣ[i]], [0,yᵣ[i]], label = "")   #рисуем две точки в координатах 0, 0 и в точке с грузом
	plot!([0,xᵣ[i]], [0,yᵣ[i]], label = "")    #рисуем нить, соединяющую точки
	plot!([xᵣ[i], xᵣ[i]+ υₓ[i] * 0.1], [yᵣ[i], yᵣ[i] + υᵧ[i]*0.1], label = "", line = (:arrow, 0.5, 2, :red))    #Для красоты нарисуем вектор скорости, масштаб скорости 0.1 относительно координат
	#
	plot!(xlims=(-L-0.1,L+ 0.1), ylims = (-L -0.1, L + 0.1), aspect_ratio =1)  #задаем диапазон для рисования графика от -L  до L и равенство масштабов
	title!("τ = $(τᵣ[i])")    #в заголовке пишем текущее время
end

Animation("/tmp/jl_uLPXik", ["000001.png", "000002.png", "000003.png", "000004.png", "000005.png", "000006.png", "000007.png", "000008.png", "000009.png", "000010.png"  …  "000292.png", "000293.png", "000294.png", "000295.png", "000296.png", "000297.png", "000298.png", "000299.png", "000300.png", "000301.png"])

gif(anim, "anim_fps15.gif", fps = 15)   #соберем все рисунки в одну анимацию

[ Info: Saved animation to /user/anim_fps15.gif

Методы решения инженерных задач

Решение дифференциальных уравнений

Численное решение дифференциальных уравнений

Решение систем уравнений

Дифференциальные уравнения высших порядков

Краевая задача

Пример решения задачи из механики

Рассмотрим пример решения задачи маятника.

Тип

Краткое описание

Связанные материалы

Теги

Категории

Языки

Форматы

Уровень

Статус

Опубликовано

Обновлено

Источник