Если у вас возникли сложности с курсовой, контрольной, дипломной, рефератом, отчетом по практике, научно-исследовательской и любой другой работой - мы готовы помочь.
Предоплата всего
Подписываем
Содержание
1. Постановка задачи
. Исследование методов решения
. Аналитическое решение задачи
. Численное решение
.1 Исходные данные
.2 Описание функций программы
.3 Исходный код
.4 Отладка и тестирование программы
.5 Исследование
. Выводы
Построить таблицу и график решения линейного дифференциального уравнения, заданным численным методом и аналитически, а также с использованием стандартной процедуры ode45.
Дифференциальное уравнение:
0 y’’’(t) + a1 y’’(t) + a2 y’(t) + y(t) = 10.
Начальные условия:
(0) = y’(0) = y’’(0) = 0.
Интервал:
Є [0, tm].
Сравнить полученные результаты с аналитическим решением. Исследовать зависимость погрешности решения от выбора шага интегрирования.
Исходные данные:
0 = 2;1 = 1.02;2 = 2.01;m = 20 (с).
В ходе работы будет использоваться метод Адамса-Башфорта 3-го порядка.
Общая формула Адамса-Башфорта:
Метод Адамса-Башфорта 3-го порядка является трехшаговым (поскольку использует информацию в трех точках), явным и имеет формулу:
.
Общая форма ЛДУ:
ЛДУ уравнение можно привести к системе, а именно:
. пусть ,
где - неизвестные числа;
2. подставим выше указанную замену в ЛДУ и выразим оттуда . Получим:
;
3. в результате получаем систему ЛДУ(СЛДУ):
Полученную СЛДУ можно переписать в матричном виде:
, где ; ;
= f, следовательно можно записать для конкретного заданного уравнения:
; ; ; n = 3
количество узлов.
Для расчёта дополнительных начальных условий требуется стартовый метод, например, одношаговый метод Рунге-Кутты.
k+1 = xk+h*(k1+2*k2+2*k3+k4)/6, где= f(tk, xk)= f(tk+h/2, xk+h*k1/2)= f(tk+h/2, xk+h*k2/2)= f(tk+h, xk+h*k3)= (tk -t0 )/n
y’’’(t) + 1.02 y’’(t) + 2.01 y’(t) + y(t) = 10
Из заданного дифференциального уравнения найдем корни характеристического многочлена.
С помощью следующей программы:
= [0 1 0; 0 0 1; -1/2 -2.01/2 -1.02/2];= eig(A)
Результат:
.5000
.0050 + 1.0000i
.0050 - 1.0000i
Так как собственные значения различаются больше, чем на порядок, то из этого следует, что эта система достаточно жесткая. Кроме того, ввиду того, что вещественные части всех собственных значений отрицательны, система асимптотически устойчива.
Соответствующие им частные решения ДУ:
Общее решение нашего ДУ:
Начальные условия:
(0) = y’(0) = y’’(0) = 0.’’’(0) = 10 / a0 = 5
.m
= vpa(dsolve('2 * D3y + 1.02 * D2y + 2.01 * Dy + y = 10', 'y(0) = 0', 'D2y(0) = 0', 'Dy(0) = 0', 'x'))
Результат:= 10.0 - (1.9678714859437751004016064257028 * cos(0.99998749992187402342224094390368*x)) / exp(0.005*x) - (4.0259539391966028659938403448691 * sin(0.99998749992187402342224094390368*x)) / exp(0.005*x) - 8.0321285140562248995983935742972 / exp(0.5*x)
) Дано:
0 = 2;1 = 1.02;2 = 2.01;m = 20 (с);= [0;0;0].
) Приведение ДУ к форме Коши
) Расчет общего решения ДУ Y(k) аналитическим методом
) Расчет с помощью ode45
) Расчет первых трёх значений с помощью метода Рунге-Кутты
) Расчет с помощью метода Адамса-Башфорта 3 -го порядка
) Вывод графиков расчета разными методами на экран
дифференциальный уравнение интегрирование программа
Основные функции и переменные, использованные в реализованной программе:kursach2(a, y0, t0, tm, h) - функция, параметрами которой являются: коэффициенты в уравнении, начальные условия ЛДУ, интервал времени, шаг= R(t,a)- функция расчета методом Рунге-Кутты
[T,Y] = AB(tspan,a)- - функция, осуществляющая стартовый метод; выходными параметрами являются вектор-столбец Y и параметр T; входными параметрами являются начало и конец отрезка интегрирования; исходная матрица а
kursach2(a, y0, t0, tm, h)= [2, 1.02, 2.01]; % коэффициенты ДУ= [ 0, 1, 0; 0, 0, 1; - 1/a(1), - a(3)/a(1), - a(2)/a(1)]; % матрица A= 0;= 20;= [0; 0; 10/a(1)];= 1; % шаг= [0; 0; 0]; % начальные условия= (tm - t0)/h; % количество точекk = 1:n + 1 % цикл по всем точкам для задания аналитического решения(k) = t0 + (k - 1) * h;(k) = 10. -8.03212851405622489959839357430 * exp(-.5 * T(k)) - 4.02595393919660286599384034488 * exp(-.5e-2 * T(k)) * sin(.999987499921874023422240943905 * T(k)) - 1.96787148594377510040160642570 * exp(- .500000000000000000000000000000e-2 * T(k)) * cos(.999987499921874023422240943905 * T(k));= 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
% построение графиков(1);(2, 1, 1);(T, (Y1(:, 1)));on;('Решение ode45');(2, 1, 2);(T, (Y' - Y2(:, 1)));on;('Погрешность ode45');(2);(2, 1, 1);(T, (Y' - Y1(:, 1)));on;('Метод Адамсa-Башвортa');(2, 1, 2);(T, (Y2(:, 1)));on;('Погрешность Адамса-Башворта');(3);(T, Y);on;('Аналитическое');X1 = R(t, a) % Рунге= [ 0, 1, 0; 0, 0, 1; - 1/a(1), - a(3)/a(1), - a(2)/a(1)];= [0; 0; 10/a(1)];= [0; 0; 0]; % начальные условия= A * y + f;(:, 1) = y;= t; %k = 1:2 % цикл= A * y + f; %= A * (y + h/2 * k1) + f;= A * (y + h/2 * k2) + f;= A * (y + h * k3) + f;= y + h * (k1 + 2 * k2 + 2 * k3 + k4)/6;= A * y + f;(:, k + 1) = y;% конец
% явный метод адамса - башфорта[T, Y] = AB(tspan, a) % диапазон времени= length(tspan); % кол - во шагов и точек разбиения= tspan(n)/(n - 1); % шаг
% X(:, 1) = y0'; % начальные значения_own = (R(h, a)); % рассчитываем с помощью метода Рунге - Кутты первые три значения(:, 1) = f_own(:, 1);(:, 2) = f_own(:, 2);(:, 3) = f_own(:, 3); % записываем их в столбцыk = 3:n - 1(:, 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)) ); % метод адамса - башфорта= tspan; Y = X';
Для сравнения погрешностей при разных шагах и анализа устойчивости, будут проводиться вычисления с различным шагом.
Найдём h критическое для метода Адамса-Башфорта 3го порядка:*| л | < 6/11
л = -0.5< 6/(11*| л |)
критическое ~1,1критическое получилось больше единицы, поэтому критическим шагом в данной работе будем считать h = 1.
Ниже представлены графики для различных h и анализ погрешностей.= 0,005
= 0,01
h = 0,05
= 0,1
= 0,5
= 1
|
Шаг h |
Погрешность метода А-Б |
Погрешность ode45 |
|
0,005 |
3*10-6 |
2*10-5 |
|
0,01 |
2,5*10-5 |
2*10-5 |
|
0,05 |
3*10-3 |
2*10-5 |
|
0,1 |
0,03 |
2*10-5 |
|
0,5 |
2,5 |
2*10-5 |
|
1 |
1,3*104 |
2*10-5 |
После проведения тестирования видно, что погрешность метода Адамса-Башфорта увеличивается с увеличением ширины шага. На шаге h < 0,01 она становится меньше погрешности стандартного ode45.
Можно сделать вывод, что при равномерном выборе величины шага метод Адамса-Башфорта является недостаточно устойчивым для большой величины шага.
Точности решения можно добиться, если подбирать величину шага с течением времени, обращая внимание на требования устойчивости.
В ходе работы была реализована программа решения ДУ методом Адамса-Башфорта 3-го порядка, использующая также метод Рунге-Кутты.
Была проведена отладка и тестирование программы, где полученное с помощью данного метода решение было сравнено с решением с помощью встроенного ode45. В результате работы была найдена зависимость погрешности интегрирования от выбранного шага.