4.3. Приближенное решение обыкновенных дифференциальных уравнений первого порядка

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

Обыкновенным дифференциальным уравнением (ОДУ) n-го порядка называется следующее уравнение, которое содержит одну или несколько производных от искомой функции y(x):

где обозначает производную порядка n некоторой функции y(x), x – это независимая переменная.

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

Общее решение ОДУ n-го порядка содержит n произвольных констант C1, C2, …, Cn

Частное решение ОДУ получается из общего, если константам интегрирования придать некоторые значения, определив некоторые дополнительные условия, количество которых позволяет вычислить все неопределенные константы интегрирования.

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

Численное решение ДУ (частное) заключается в вычислении функции y(x) и ее производных в некоторых заданных точках x₁,x₂,…,x_N, лежащих на определенном отрезке.

4.3.1. Метод Эйлера

Простейшим численным методом решения обыкновенных дифференциальных уравнений является метод Эйлера. В его основе лежит аппроксимация производной отношением конечных приращений зависимой (y) и независимой (x) переменных между узлами равномерной сетки:

,

где y_i+1 – это искомое значение функции в точке x_i+1.

Если теперь преобразовать это уравнение, и учесть равномерность сетки интегрирования, то получится итерационная формула, по которой можно вычислить y_i+1, если известно y_i в точке х_i:

Рис.4.1. Графическая иллюстрация метода Эйлера

Таким образом, суть метода Эйлера заключается в замене функции y(x) на отрезке интегрирования прямой линией, касательной к графику в точке x=xi. Если искомая функция сильно отличается от линейной на отрезке интегрирования, то погрешность вычисления будет значительной. Ошибка метода Эйлера прямо пропорциональна шагу интегрирования:

Пример 4.3.1 Рассчитать кинетику гомогенной химической реакции, используя метод Эйлера.

VAR ca0, cb0, cc0, cd0, k1, k2, k3, k4, t0, tk, h, ca, cb, cc, cd: Real;

{Описание переменных концентрации, констант скорости, время реагирования, шага}

f1, f2 :Text;

{Файловые переменные}

BEGIN

assign(f1,'al_ar.in');

{Файл ввода данных}

assign(f2,'al_ar.out');

{Файл вывода расчетов}

reset(f1);

readln(f1,ca0,cb0,cc0,cd0,k1,k2,k3,k4,t0,tk,h);

{ Ca0, Cb0, Cc0, Cd0 - Концентрации веществ для t0 }

{ k1, k2, k3, k4 - Константы скорости для каждой реакции }

{ t0, tk - Начальное и конечное время реагирования }

{ Ca, Cb, Cc, Cd - Переменные для расчета концентраций }

close(f1);

ca:=ca0; cb:=cb0; cc:=cc0; cd:=cd0;

rewrite(f2);

writeln(f2,'Метод Эйлера');

writeln(f2,'Время Концентрация');

writeln(f2,' А B С D ');

REPEAT

ca:=ca0+h*((-1*(k1+k4))*ca0+k3*cd0);

cb:=cb0+h*(k1*ca0-k2*cb0);

cc:=cc0+h*k2*cb0;

cd:=cd0+h*((k4*ca0)-(k3*cb0));

writeln(f2,' ',t0:8:4,' ',ca:8:4,' ',cb:8:4,' ',cc:8:4,' ',cd:8:4,' ');

ca0:=ca; cb0:=cb; cc0:=cc; cd0:=cd; t0:=t0+h;

UNTIL t0>=tk;

writeln(f2,' ',t0:8:4,' ',ca:8:4,' ',cb:8:4,' ',cc:8:4,' ',cd:8:4,' ');

close(f2);

END.