Математический анализ
Численное решение обыкновенных дифференциальных уравнений
Обыкновенные дифференциальные уравнения первого порядка
При решении различных задач математики, физики, химии и других наук часто пользуются математическими моделями в виде уравнений, связывающих независимую переменную, искомую функцию и ее производные. Такие уравнения называются дифференциальными уравнениями (ДУ).
Если искомая функция зависит от одной независимой переменной , то ДУ называют обыкновенным дифференциальным уравнением (ОДУ). Если искомая функция зависит от нескольких переменных, то ДУ называют дифференциальным уравнением в частных производных. В этом курсе мы будем рассматривать только ОДУ.
Наивысший порядок производной искомой функции, входящей в ДУ, называется порядком этого уравнения.
ОДУ первого порядка в общем случае можно записать в виде
Если из ОДУ можно выразить производную от искомой функции , то его называют ОДУ первого порядка, разрешенным относительно производной:
Общее и частное решение дифференциального уравнения
Решением ОДУ называется функция , которая при подстановке в уравнение обращает его в тождество.
Любое ДУ имеет бесконечное множество решений, отличающихся друг от друга постоянными величинами. Например, решением уравнения является любая функция вида , где – произвольная константа.
Общим решением ОДУ первого порядка называется функция , содержащая одну произвольную постоянную. Эта функция является решением ДУ при любом фиксированном значении .
Частным решением ОДУ первого порядка называется любая функция , полученная из общего решения при конкретном значении постоянной .
Чтобы найти какое-либо частное решение ДУ, на общее решение необходимо наложить
начальное условие: при заданном значении независимой переменной функция должна быть равна заданному числу :
Задача отыскания частного решения ОДУ первого порядка, удовлетворяющего заданному начальному условию, называется задачей Коши.
Пример ОДУ первого порядка: Модель движения ракеты
Найдем зависимость высоты взлетающей ракеты над уровнем земли от времени . Пусть ракета взлетает с постоянным ускорением , где – ускорение свободного падения. Тогда скорость ракеты, с одной стороны, равна произведению ускорения на время, т.е. , а с другой стороны, скорость равна производной от координаты (высоты) по времени . Приравнивая правые части этих выражений, получим ОДУ первого порядка относительно неизвестной функции :
Начальное условие имеет следующий вид: ракета находится на нулевой высоте () в начальный момент времени (), т.е.
Данное уравнение легко решается аналитически. Для этого перенесем в правую часть уравнения и проинтегрируем обе части:
Отсюда получаем общее решение:
Подставим начальное условие в общее решение: . Отсюда
Подставляя значение константы в общее решение, получим искомое частное решение:
Численные методы решения ОДУ
Методами математического анализа для ряда ОДУ можно получить точное (аналитическое) решение в виде формулы . Однако для многих ОДУ аналитическое решение оказывается слишком сложным или его вообще невозможно получить. Например, ОДУ не решается аналитически. В этих случаях используют приближенные (численные) методы решения ОДУ. Численное решение задачи Коши представляет собой таблицу значений независимой переменной , называемых узлами, и соответствующих им значений искомой функции .
| ... | |||
|---|---|---|---|
| ... |
Возьмем равноотстоящие узлы: . Приближенные значения искомой функции в узлах можно найти несколькими численными методами. Приведем простейшие из них.
Метод Эйлера:
Метод Эйлера–Коши:
Метод Рунге–Кутта:
где
Метод Эйлера имеет первый порядок точности, метод Эйлера–Коши – второй порядок точности, метод Рунге–Кутта – четвертый порядок точности.
Непосредственная реализация методов численного решения ОДУ в Engee
Решим задачу Коши , на отрезке методами Эйлера и Эйлера–Коши. Разобьем отрезок на частей.
using Plots;
u0 = 1;
a = 0;
b = 1;
n = 100;
h = (b - a) / n;
t = zeros(n+1);
func1 = zeros(n+1);
func2 = zeros(n+1);
u1 = zeros(n+1);
u2 = zeros(n+1);
t[1] = a;
u1[1] = u0;
u2[1] = u0;
for i in 1:n
t[i+1] = t[i] + h;
func1[i] = u1[i] * t[i];
u1[i+1] = u1[i] + h * func1[i]; # метод Эйлера
func2[i] = u1[i+1] * t[i];
u2[i+1] = u2[i] + h/2 * (func2[i] + func2[i]) # метод Эйлера--Коши
end
Построим графики приближенных решений, полученных двумя численными методами, и график точного решения .
u_acc = exp.((t.^2) / 2)
plot(t, u1, label="методом Эйлера");
plot!(t, u2, label="методом Эйлера--Коши")
plot!(t, u_acc, label="точное решение")
✏️Задание 1
Найдите в Engee решение ОДУ с начальным условием на отрезке методами Эйлера и Эйлера–Коши. Постройте графики приближенных решений и сравние их с графиком точного решения .
Решение
using Plots;
u0 = 1;
a = 0;
b = 2;
n = 100;
h = (b - a) / n;
t = zeros(n+1);
func1 = zeros(n+1);
func2 = zeros(n+1);
u1 = zeros(n+1);
u2 = zeros(n+1);
t[1] = a;
u1[1] = u0;
u2[1] = u0;
for i in 1:n
t[i+1] = t[i] + h;
func1[i] = u1[i] + t[i];
u1[i+1] = u1[i] + h * func1[i]; # метод Эйлера
func2[i] = u1[i+1] + t[i];
u2[i+1] = u2[i] + h/2 * (func2[i] + func2[i]) # метод Эйлера--Коши
end
u_acc = 2*exp.(t) .- t .- 1; # точное решение
# построение графиков
plot(t, u1, label="методом Эйлера");
plot!(t, u2, label="методом Эйлера--Коши")
plot!(t, u_acc, label="точное решение")
Численное решение ОДУ с помощью встроенных функций Engee
Пример. Найдем в Engee частное решение ОДУ первого порядка с начальным условием на отрезке .
Для этого воспользуемся функциями ODEProblem и solve, находящимися в библиотеке .DifferentialEquations. Сначала зададим правую часть уравнения f (предполагается, что в левой части уравнения находится производная от искомой функции u'), начальное значение искомой функции u0 и отрезок независимой переменной tspan, на котором мы будем искать решение.
using Plots, DifferentialEquations;
f(u,p,t) = 0.98*u;
u0 = 1.0;
tspan = (0.0, 1.0);
Затем с помощью функции ODEProblem составим задачу Коши. Аргументами этой функции являются правая часть дифференциального уравнения f, начальное значение искомой функции u0 и отрезок tspan.
prob = ODEProblem(f, u0, tspan);
И, наконец, найдем частное решение задачи Коши с помощью функции solve.
sol = solve(prob)
Построим график найденного решения.
plot(sol)
✏️Задание 2
Найдите в Engee решение ОДУ с начальным условием на отрезке с помощью встроенных функций Engee. Постройте график приближенного решения и сравние его с графиком точного решения .
Решение
using Plots, DifferentialEquations;
f(u,p,t) = u + t;
u0 = 1.0;
tspan = (0.0, 2.0);
prob = ODEProblem(f, u0, tspan);
sol = solve(prob)
plot(sol, label="приближенное решение")
t = collect(0:(2-0)/100:2);
u_acc = 2*exp.(t) .- t .- 1;
plot!(t, u_acc, label="точное решение")