Методы решения обыкновенных дифференциальных уравнений

Методы решения обыкновенных дифференциальных уравнений

Решить дифференциальное уравнение (1) для вертикальных колебаний под действием вынуждающей силы.

(1).

Где x — отклонение системы от исходного положения, t — время, m — масса, в — коэффициент вязкого трения, k — коэффициент жесткости амортизаторов, щ — циклическая частота.

Решим это дифференциальное уравнение для следующих параметров:

x = 0, dx/dt = 0, t = 0, в=0.5кг/с, k = 5 Н/м, m = 1 кг, Fm = 2000Н,.

щ = 0.3 рад/с.

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

Метод Эйлера является явным, одношаговым методом первого порядка точности, основанном на аппроксимации интегральной кривой кусочно-линейной функцией, т. н. ломаной Эйлера. В общем случае сводиться к уравнению (2), где h — шаг. Для уравнения 2 порядка необходимо дважды вычислить данное уравнение.

(2).

Метод легко реализуем, и требует минимум вычислений, но в теории дает низкую точность.

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

Приближенное значение в точках:

Коэффициенты:

Сначала решим дифференциальное уравнение с помощью функции ode45(Метод Рунге-Кутты 4, 5 порядка). Создадим файл-функцию для нашего уравнения.

functiondx = diff_fun (t, y).

% Константы.

m=1;

B=0.5;

k=5;

w=0.3;

%Уравнения.

dx=[y (2); 1/m*(2000*abs (cos (w*t))-B*y (2)-k*y (1))];

Запишем скрипт для нахождения решения.

clc; clear; closeall;

Y0=[0 0]; % Начальные значения.

h=0.01; % Шаг.

t0=0;

t1=30;

T0=[t0:h:t1]; % Интервал для решения.

[T, Y] = ode45('diff_fun', T0, Y0);

holdon.

plot (T, Y (, 1), 'b-', 'LineWidth', 2).

plot (T, Y (, 2), 'r-', 'LineWidth', 2).

gridon.

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

clc; clear; closeall;

Y0=[0 0]; % Начальные значения.

h=0.01; % Шаг.

t0=0;

t1=30;

%% Встроенная функция (Метод Рунге-Кутты).

T0=[t0:h:t1]; % Интервал для решения.

[T, Y] = ode45('diff_fun', T0, Y0);

hold on.

plot (T, Y (, 1), 'b-', 'LineWidth', 2).

plot (T, Y (, 2), 'r-', 'LineWidth', 2).

grid on.

Y11=Y (, 1); Y22=Y (, 2);

clearT0TY% Очищаем переменные.

%% Собственная функция (Метод Рунге-Кутты).

Y1=[]; Y2=[]; % Задаем массив.

y1=Y0(1); y2=Y0(2); % Задаем начальные значения.

T=t0:h:t1; % Задаем интервал.

for t = t0: h:t1 % Решаем.

k1=h*dx (t, y1, y2, '1');

k2=h*dx (t+h/2, y1+k½, y2+k½, '1');

k3=h*dx (t+h/2, y1+k2/2, y2+k2/2, '1');

k4=h*dx (t+h, y1+k3, y2+k3, '1');

y1=y1+(k1+2*k2+2*k3+k4)/6;

Y1=[Y1 y1];

k1=h*dx (t, y1, y2, '2');

k2=h*dx (t+h/2, y1+k½, y2+k½, '2');

k3=h*dx (t+h/2, y1+k2/2, y2+k2/2, '2');

k4=h*dx (t+h, y1+k3, y2+k3, '2');

y2=y2+(k1+2*k2+2*k3+k4)/6;

Y2=[Y2 y2];

end.

plot (T, Y1, 'm-', 'lineWidth', 2).

plot (T, Y2, 'g-', 'LineWidth', 2).

legend ('y1 — ode45', 'y2 — ode45', 'y1 — МетодРунге-Кутты', 'y2 — МетодРунге-Кутты').

xlabel ('Время, с');

ylabel ('Амплитуда');

clearY2Y1.

%% Собственная функция (Метод Эйлера).

Y1=[]; Y2=[]; % Задаем массив.

y1=Y0(1); y2=Y0(2); % Задаем начальные значения.

for t = t0: h:t1.

y1=y1+h*dx (t, y1, y2, '1');

y2=y2+h*dx (t, y1, y2, '2');

Y1=[Y1 y1];

Y2=[Y2 y2];

end.

figure.

hold on; grid on.

plot (T, Y11, 'b-', 'LineWidth', 2).

plot (T, Y22, 'r-', 'LineWidth', 2).

plot (T, Y1, 'm-', 'lineWidth', 2).

plot (T, Y2, 'g-', 'LineWidth', 2).

legend ('y1 — ode45', 'y2 — ode45', 'y1 — МетодЭйлера', 'y2 — МетодЭйлера').

xlabel ('Время, с');

ylabel ('Амплитуда');

Функция:

function y = dx (t, y1, y2, n).

% Константы.

m=1;

B=0.5;

k=5;

w=0.3;

switch n.

case {'1'}.

y=y2;

case {'2'}.

y=1/m*(2000*abs (cos (w*t))-B*y2-k*y1);

end.

Получаем графики на рисунке 1, 2. Видим, что решение по методу Рунге-Кутты 4 порядка незначительно отличается от решения через функцию ode45. Но совершенно неожиданные результаты дает метод Эйлера.

(см. рисунок 2), данные полученные этим методом практически идеально совпадают с данными функции ode45.

Рисунок 1 — Сравнение функции ode45 и метода Рунге-Кутты 4 порядка.

Рисунок 2 — Сравнение ode45 и метода Эйлера.

Заключение

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