Интегрирование обыкновенных дифференциальных уравнений

• 1. Постановка задачи
• 2. Исследование методов решения
• 3. Аналитическое решение задачи
• 4. Численное решение
• 4.1 Исходные данные
• 4.2 Описание функций программы
• 4.3 Исходный код
• 4.4 Отладка и тестирование программы
• 4.5 Исследование
• 5. Выводы

1. Постановка задачи

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

Дифференциальное уравнение:

a₀ y'''(t) + a₁ y''(t) + a₂ y'(t) + y (t) = 10.

Начальные условия:

y (0) = y'(0) = y''(0) = 0.

Интервал:

t Є [0, t_m].

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

Исходные данные:

a₀ = 2;

a₁ = 1.02;

a₂ = 2.01;

t_m = 20 ©.

В ходе работы будет использоваться метод Адамса-Башфорта 3-го порядка.

2. Исследование методов решения

Общая формула Адамса-Башфорта:

Метод Адамса-Башфорта 3-го порядка является трехшаговым (поскольку использует информацию в трех точках), явным и имеет формулу:

.

Общая форма ЛДУ:

ЛДУ уравнение можно привести к системе, а именно:

1. пусть ,

где — неизвестные числа;

2. подставим выше указанную замену в ЛДУ и выразим оттуда. Получим:

;

3. в результате получаем систему ЛДУ (СЛДУ):

Полученную СЛДУ можно переписать в матричном виде:

где; ;

Bu = f, следовательно можно записать для конкретного заданного уравнения:

;;; n = 3

количество узлов.

Для расчёта дополнительных начальных условий требуется стартовый метод, например, одношаговый метод Рунге-Кутты.

x_k+1 = x_k+h*(k1+2*k2+2*k3+k4)/6, где

k1 = f (t_k, x_k)

k2 = f (t_k+h/2, x_k+h*k½)

k3 = f (t_k+h/2, x_k+h*k2/2)

k4 = f (t_k+h, x_k+h*k3)

h = (t_k -t₀ )/n

3. Аналитическое решение задачи

2y'''(t) + 1.02 y''(t) + 2.01 y'(t) + y (t) = 10

Из заданного дифференциального уравнения найдем корни характеристического многочлена.

С помощью следующей программы:

A = [0 1 0; 0 0 1; -½ -2.01/2 -1.02/2];

ei1 = eig (A)

Результат:

— 0.5000

— 0.0050 + 1.0000i

— 0.0050 — 1.0000i

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

Соответствующие им частные решения ДУ:

Общее решение нашего ДУ:

Начальные условия:

y (0) = y'(0) = y''(0) = 0.

y'''(0) = 10 / a₀ = 5

F.m

y = vpa (dsolve ('2 * D3y + 1.02 * D2y + 2.01 * Dy + y = 10', 'y (0) = 0', 'D2y (0) = 0', 'Dy (0) = 0', 'x'))

Результат:

y = 10.0 — (1.9 678 714 859 437 750 547 611 683 651 584 * cos (0.99 998 749 992 187 396 221 929 246 425 088*x)) / exp (0.005*x) — (4.259 539 391 966 028 665 293 077 741 568 * sin (0.99 998 749 992 187 396 221 929 246 425 088*x)) / exp (0.005*x) — 8.321 285 140 562 249 017 915 867 136 000 / exp (0.5*x)

4. Численное решение

4.1 Исходные данные

1) Дано:

a₀ = 2;

a₁ = 1.02;

a₂ = 2.01;

t_m = 20 (с);

y0 = [0;0;0].

2) Приведение ДУ к форме Коши

3) Расчет общего решения ДУ Y (k) аналитическим методом

4) Расчет с помощью ode45

5) Расчет первых трёх значений с помощью метода Рунге-Кутты

6) Расчет с помощью метода Адамса-Башфорта 3 -го порядка

7) Вывод графиков расчета разными методами на экран

4.2 Описание функций программы

дифференциальный уравнение интегрирование программа

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

function kursach2(a, y0, t0, tm, h) — функция, параметрами которой являются: коэффициенты в уравнении, начальные условия ЛДУ, интервал времени, шаг

X1 = R (t, a) — функция расчета методом Рунге-Кутты

[T, Y] = AB (tspan, a) — - функция, осуществляющая стартовый метод; выходными параметрами являются вектор-столбец Y и параметр T; входными параметрами являются начало и конец отрезка интегрирования; исходная матрица а

4.3 Исходный код

function kursach2(a, y0, t0, tm, h)

a = [2, 1.02, 2.01]; % коэффициенты ДУ

A = [ 0, 1, 0; 0, 0, 1; - 1/a (1), — a (3)/a (1), — a (2)/a (1)]; % матрица A

t0 = 0;

tm = 20;

f = [0; 0; 10/a (1)];

h = 1; % шаг

y0 = [0; 0; 0]; % начальные условия

n = (tm — t0)/h; % количество точек

for k = 1: n + 1% цикл по всем точкам для задания аналитического решения

T (k) = t0 + (k — 1) * h;

Y (k) = 10. -8.3 212 851 405 622 490 179 158 671 360 * exp (-.5 * T (k)) — 4.2 595 393 919 660 286 455 018 684 416 * exp (-.5e-2 * T (k)) * sin (.999 987 499 921 873 967 848 791 998 464 * T (k)) — 1.96 787 148 594 377 501 253 992 185 856 * exp (- .49 999 999 999 999 995 716 575 428 6080e-2 * T (k)) * cos (.999 987 499 921 873 967 848 791 998 464 * T (k));

end

options = odeset ('RelTol', 1e-10, 'AbsTol', 1e-5); % точность ode45

[T1, Y1] = ode45('F', T, y0, options); % решение ode45

%

[T2, Y2] = AB (T, a); % метод Адамса — Башфорта для шага h

[T3, Y3] = AB (T/2, a); % метод Адамса — Башфорта для шага h/2

% построение графиков

figure (1);

subplot (2, 1, 1);

plot (T, (Y1(, 1)));

grid on;

title ('Решение ode45');

subplot (2, 1, 2);

plot (T, (Y' - Y2(, 1)));

grid on;

title ('Погрешность ode45');

figure (2);

subplot (2, 1, 1);

plot (T, (Y' - Y1(, 1)));

grid on;

title ('Метод Адамсa-Башвортa');

subplot (2, 1, 2);

plot (T, (Y2(, 1)));

grid on;

title ('Погрешность Адамса-Башворта');

figure (3);

plot (T, Y);

grid on;

title ('Аналитическое');

function X1 = R (t, a) % Рунге

A = [ 0, 1, 0; 0, 0, 1; - 1/a (1), — a (3)/a (1), — a (2)/a (1)];

f = [0; 0; 10/a (1)];

y = [0; 0; 0]; % начальные условия

y1 = A * y + f;

X1(, 1) = y;

h = t; %

for k = 1:2% цикл

k1 = A * y + f; %

k2 = A * (y + h/2 * k1) + f;

k3 = A * (y + h/2 * k2) + f;

k4 = A * (y + h * k3) + f;

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

y1 = A * y + f;

X1(, k + 1) = y;

end % конец

% явный метод адамса — башфорта

function [T, Y] = AB (tspan, a) % диапазон времени

n = length (tspan); % кол — во шагов и точек разбиения

h = tspan (n)/(n — 1); % шаг

% X (, 1) = y0'; % начальные значения

f_own = (R (h, a)); % рассчитываем с помощью метода Рунге — Кутты первые три значения

X (, 1) = f_own (, 1);

X (, 2) = f_own (, 2);

X (, 3) = f_own (, 3); % записываем их в столбцы

for k = 3: n — 1

X (, k + 1) = X (, k) + h/12 * (23 * F (tspan (k), X (, k)) — 16 * F (tspan (k — 1), X (, k — 1)) + 5 * F (tspan (k — 2), X (, k — 2))); % метод адамса — башфорта

end

T = tspan; Y = X';

4.4 Отладка и тестирование программы

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

Найдём h критическое для метода Адамса-Башфорта 3го порядка:

h*| л | < 6/11

л = -0.5

h < 6/(11*| л |)

h критическое ~1,1

h критическое получилось больше единицы, поэтому критическим шагом в данной работе будем считать h = 1.

Ниже представлены графики для различных h и анализ погрешностей.

h = 0,005

h = 0,01

h = 0,05

h = 0,1

h = 0,5

h = 1

4.5 Исследование

После проведения тестирования видно, что погрешность метода Адамса-Башфорта увеличивается с увеличением ширины шага. На шаге h < 0,01 она становится меньше погрешности стандартного ode45.

Можно сделать вывод, что при равномерном выборе величины шага метод Адамса-Башфорта является недостаточно устойчивым для большой величины шага.

Точности решения можно добиться, если подбирать величину шага с течением времени, обращая внимание на требования устойчивости.

5. Выводы

В ходе работы была реализована программа решения ДУ методом Адамса-Башфорта 3-го порядка, использующая также метод Рунге-Кутты.

Была проведена отладка и тестирование программы, где полученное с помощью данного метода решение было сравнено с решением с помощью встроенного ode45. В результате работы была найдена зависимость погрешности интегрирования от выбранного шага.