Математический анализ
Численное решение систем обыкновенных дифференциальных уравнений
Системы обыкновенных дифференциальных уравнений первого порядка
Для решения многих задач математики, физики, химии, техники нередко требуется несколько функций. Нахождение этих функций может привести к нескольким ДУ, образующим систему.
Системой ДУ называется совокупность ДУ, каждое из которых содержит независимую переменную, искомые функции и их производные.
Рассмотрим систему ДУ первого порядка, разрешенных относительно производной:
Такая система называется нормальной системой ДУ. При этом число уравнений равно числу искомых функций.
Решением системы ДУ называется совокупность из функций , удовлетворяющих каждому из уравнений этой системы.
Общее решение такой системы ДУ содержит произвольных постоянных:
Решение, получающееся из общего при конкретных значениях постоянных , называется частным решением системы.
Начальные условия для этой системы имеют вид:
Задача Коши ставится следующим образом: найти частное решение системы ДУ, удовлетворяющее заданным начальным условиям.
Численные методы решения систем ОДУ первого порядка
Для многих систем ОДУ получить аналитическое решение оказывается слишком сложно или вообще невозможно. В таких случаях применяют численные методы решения систем ОДУ. Рассмотренные в предыдущем разделе методы решения ОДУ первого порядка можно использовать и для систем ОДУ первого порядка. При этом расчетные формулы претерпевают минимальные изменения: следует лишь заменить числа , функцию и начальное условие на соответствующие векторы , и .
Например, расчетная формула метода Эйлера примет вид:
Покоординатная запись этого выражения выглядит так:
Здесь первый индекс обозначает номер искомой функции, а второй индекс – номер узла, в котором вычисляется значение функции.
Аналогично можно записать формулы для решения систем ОДУ методами Эйлера–Коши и Рунге–Кутта.
Задача Коши для систем ОДУ является жесткой, т.е. несмотря на медленное изменение искомых функций, расчет приходится вести с очень мелким шагом. Увеличение шага приводит к катастрофически большому росту погрешности.
Непосредственная реализация методов численного решения систем ОДУ в Engee
Рассмотрим систему из двух ОДУ первого порядка с двумя неизвестными функциями и :
Начальные условия имеют вид: , .
Запишем алгоритм метода Эйлера для численного решения такой системы:
Пример. Решим в Engee методом Эйлера систему ОДУ
при начальных условиях , на отрезке .
Разобьем отрезок на частей и реализуем алгоритм метода Эйлера аналогично тому, как мы это делали в предыдущем разделе для одиночных ОДУ. Затем построим графики решения.
using Plots;
a = 0;
b = 5;
n = 1000;
h = (b - a) / n;
t = zeros(n+1);
func1 = zeros(n+1);
func2 = zeros(n+1);
x = zeros(n+1);
y = zeros(n+1);
t[1] = a;
x[1] = 1;
y[1] = 2;
for i in 1:n
t[i+1] = t[i] + h;
func1[i] = 4*x[i] - 3*y[i];
x[i+1] = x[i] + h * func1[i];
func2[i] = 2*x[i] - 3*y[i];
y[i+1] = y[i] + h * func2[i];
end
plot(t, x, label="x(t)")
plot!(t, y, label="y(t)")
✏️Задание 1
Решите в Engee методом Эйлера систему уравнений
при начальных условиях , на отрезке , разбивая отрезок на частей. Постройте графики решения.
Решение
using Plots;
a = 0;
b = 1;
n = 1000;
h = (b - a) / n;
t = zeros(n+1);
func1 = zeros(n+1);
func2 = zeros(n+1);
x = zeros(n+1);
y = zeros(n+1);
t[1] = a;
x[1] = 2;
y[1] = 3;
for i in 1:n
t[i+1] = t[i] + h;
func1[i] = x[i] - y[i];
x[i+1] = x[i] + h * func1[i];
func2[i] = -4*x[i] + y[i];
y[i+1] = y[i] + h * func2[i];
end
plot(t, x, label="x(t)")
plot!(t, y, label="y(t)")
Численное решение систем ОДУ с помощью встроенных функций Engee
Для решения систем ОДУ первого порядка используются те же встроенные функции ODEProblem и solve, которые мы использовали в предыдущем разделе для решения одиночных ОДУ. Решаемая система уравнений должна представлять собой нормальную систему. При этом задача Коши записывается в векторной форме: правые части уравнений и начальные условия представляют собой векторы.
Пример. Решим в Engee систему из трех нелинейных ОДУ, играющую большую роль в теории хаоса и называемую аттрактором Лоренца:
при следующих значениях параметров: , , и начальных условиях: , , .
Сначала зададим вектор-функцию lorenz, представляющую собой вектор правых частей уравнений системы. Здесь u – вектор неизвестных функций x, y и z; du – вектор производных от неизвестных функций dx/dt, dy/dy и dz/dt; а p – вектор параметров , , .
using Plots, DifferentialEquations;
function lorenz(du, u, p, t)
sigma, rho, beta = p
du[1] = sigma*(u[2] - u[1])
du[2] = u[1]*(rho - u[3]) - u[2]
du[3] = u[1]*u[2] - beta*u[3]
end
Затем зададим вектор начальных условий u0
u0 = [1.0, 0.0, 0.0];
вектор параметров p
p = (10, 28, 8/3);
и отрезок независимой переменной , на котором мы будем искать решение:
tspan = (0.0, 100.0);
Теперь найдем решение заданной системы ОДУ с помощью тех же встроенных функций Engee, которые мы использовали в предыдущем разделе для решения одиночных ОДУ. Задача Коши задается с помощью функции ODEProblem. Ее аргументами являются вектор-функция правых частей уравнений lorenz, вектор начальных условий u0, отрезок tspan и вектор параметров p, входящих в систему уравнений. Затем с помощью функции solve находится решение задачи Коши.
prob = ODEProblem(lorenz, u0, tspan, p);
sol = solve(prob)
И, наконец, построим графики решения:
plot(sol)
✏️Задание 2
Решите с помощью встроенных функций Engee систему уравнений
при начальных условиях , на отрезке . Постройте графики решения.
Решение
using Plots, DifferentialEquations;
function f(du, u, t, p)
du[1] = u[1] - u[2]
du[2] = -4*u[1] + u[2]
end
u0 = [2.0, 3.0];
tspan = (0.0, 1.0);
prob = ODEProblem(f, u0, tspan);
sol = solve(prob)
plot(sol)