Численное решение задачи Коши
Оглавление
Аннотация.
Оглавление.
. Задача № 1 (1.4)
.1 Постановка задачи.
.2 Исходные данные.
.3 Решение поставленной
задачи.
. Задача № 2 (2.2)
.1 Постановка задачи.
.2 Исходные данные.
.3 Решение поставленной
задачи.
. Задача № 3 (6.2)
.1 Постановка задачи.
.2 Исходные данные.
.3 Решение поставленной
задачи.
Заключение.
Список литературы.
1. Задача № 1 (1.4)
1.1 Постановка
задачи
Найти приближенное решение задачи Коши для
обыкновенного дифференциального уравнения (ОДУ) 1 порядка
(1)
и оценить погрешность решения
задачи.
Порядок решения задачи:
1. Задать исходные данные: функцию f
правой части, начальное значение 
.
. Используя функцию eyler
(см. ПРИЛОЖЕНИЕ B), найти
приближенное решение задачи Коши с шагом h=0.1
по явному методу Эйлера.
. Используя встроенную функцию rkfixed
пакета MATHCAD,
найти приближенное решение задачи Коши с шагом h=0.1
по методу Рунге-Кутты 4 порядка точности (см. ПРИЛОЖЕНИЕ B).
. Найти решение задачи Коши аналитически.
. Построить таблицы значений приближенных и
точного решений. На одном чертеже построить графики приближенных и точного
решений.
. Оценить погрешность приближенных решений двумя
способами:) по формуле 
; здесь
и
-
значения точного и приближенного решений в узлах сетки 
,
i=1,..N;)
по правилу Рунге (по правилу двойного пересчета) (см. ПРИЛОЖЕНИЕ C).
7. Выяснить, при каком значении шага h=h*
решение, полученное по методу Эйлера, будет иметь такую же погрешность (см. п.
6а), как решение, полученное с помощью метода Рунге-Кутты с шагом h=0.1.
УКАЗАНИЕ. В п. 7 рекомендуется провести серию
вычислений решения по методу Эйлера, дробя шаг h
пополам.
1.2 Исходные
данные
|
N
|
f(t,y)
|
t0
|
T
|
y0
|
|
1.4
|
|
0
|
1
|
1
|
.3 Решение поставленной задачи
1) Задача
Коши:
y’(t)=
,
t0=0, T=1, y0=1.
Исходные данные:
Начальное значение:
Концы отрезка:
Шаг сетки:
Число узлов сетки:
2) Функция,
реализующая явный метод Эйлера, возвращает вектор решения:
Входные параметры:
f - функция правой части;
y0 - начальное значение;- начальная точка
отрезка;- шаг сетки;- число узлов сетки.
3) Приближенное
решение задачи Коши с шагом h=0.1
по методу Рунге-Кутты 4 порядка точности с помощью встроенной функции rkfixed
пакета MATHCAD.
Функция rkfixed возвращает матрицу, первый
столбец которой содержит узлы сетки, а второй - приближенное решение в этих
узлах.
4) Аналитическое
решение задачи:
,
,
=
,
,
По методу вариации произвольной постоянной
заменим постоянную С на функцию C(t) и решим неоднородное уравнение:
Подставляем в исходное уравнение:
=
,
,
Решение в MathCad:
5) Решения,
полученные различными способами:
Метод Эйлера:
Метод Рунге-Кутты:
Точное решение:
Графики приближенных и точного решений:
6) Рассчитаем
погрешность полученных приближенных решений:
Погрешность метода Эйлера:
Вычисление погрешности по правилу Рунге:
Вычисление приближенных решений с шагом h/2:
Вычисление погрешностей:
Значение погрешностей:
7) Проведём
серию вычислений решения по методу Эйлера, дробя шаг h
пополам.
Первая итерация:
Вторая итерация:
Третья итерация:
И т.д.
Девятая итерация:
При значении шага hd=h/1024=0.1/1024=0,0000977
решение, полученное по методу Эйлера, будет иметь примерно такую же
погрешность, как и решение, полученное с помощью метода Рунге-Кутты с шагом h=0.1.
Погрешность решения по методу Эйлера с шагом hd=0,0000977
:
Погрешность решения по методу Рунге-Кутты с
шагом h=0.1 :
Вывод: явный метод Эйлера - это численный метод
1-го порядка точности. Метод Рунге-Кутты - это метод 4-го порядка точности. Это
означает, что при одном и том же значении шага, метод Рунге-Кутты даёт более
точное значение. Поэтому погрешности методов сильно (на несколько порядков)
отличаются. В рассмотренном выше примере с помощью метода Рунге-Кутты было
получено решение, которое совпадает с решением, полученным аналитическим путём.
2. Задача
№ 2 (2.2)
2.1 Постановка
задачи
Задача Коши для ОДУ 2 порядка
, 
описывает движение груза массы m,
подвешенного к концу пружины. Здесь x(t)
- смещение груза от положения равновесия, H
- константа, характеризующая силу сопротивления среды, k
-коэффициент упругости пружины, f(t)
- внешняя сила. Начальные условия: - смещение груза в начальный момент времени t=0,
-
скорость груза в начальный момент времени. Промоделировать движение груза на
временном отрезке [0,T]
при заданных в индивидуальном варианте трех наборах (I,
II, III)
значений параметров задачи. Для каждого набора по найденной таблице (или
графику) решения задачи определить максимальное и минимальное значения функции x(t)
и моменты времени, в которые эти значения достигаются. Предложить свой вариант
задания параметров, при которых характер колебаний груза существенно отличается
от рассмотренного ранее.
ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:
. Заменить исходную задачу эквивалентной задачей
Коши для системы ОДУ 1 порядка:
(2)
. Для каждого варианта выбора
параметров решить задачу (2) с помощью метода Рунге-Кутты 4 порядка точности с
шагом h=0.1.
. Для каждого варианта выбора
параметров построить график найденного решения. Сравнить характер движения
груза и дать интерпретацию полученного движения.
. Для каждого варианта выбора
параметров определить требуемые в задаче характеристики.
УКАЗАНИЕ. В п. 2 использовать
встроенную функцию rkfixed пакета MATHCAD (см.
ПРИЛОЖЕНИЕ B).
2.2 Исходные
данные
|
N
|
|
H
|
k
|
m
|
f(t)
|
x0
|
v0
|
T
|
|
2.2
|
I II III
|
1 1 1
|
1 1 1
|
0.5 0.5 0.5
|
tsin(t) 0 tsin(t)
|
0 0 0
|
0 -10 -50
|
20 20 20
|
2.3 Решение поставленной задачи
набор.
Исходные данные:
Шаг сетки:
Число узлов сетки:
Формирование вектора правой части системы ОДУ и
вектора начальных условий для применения встроенной функции rkfixed:
График решения:
Найдем максимальное и минимальное значения
функции x(t):
Минимальное значение достигается в момент
времени 18.4.
Максимальное значение достигается в момент времени
15.3.
Из графика видно, что при данном наборе значений
груз совершает незатухающие колебания.
2 набор.
Исходные данные:
Шаг сетки:
Число узлов сетки:
Формирование вектора правой части системы ОДУ и
вектора начальных условий для применения встроенной функции rkfixed:
График решения:
Минимальное значение функции достигается в
момент времени 0.8.
Максимальное значение функции достигается в момент
времени 3.9.
При данном наборе значений происходит затухание
колебаний - груз останавливается.
3 набор.
Исходные данные:
Шаг сетки:
Число узлов сетки:
Формирование вектора правой части системы ОДУ и
вектора начальных условий для применения встроенной функции rkfixed:
График решения:
Минимальное значение функции достигается в
момент времени 0.8.
Максимальное значение функции достигается в
момент времени 3.9.
При данном наборе значений, как и в первом
случае, груз совершает незатухающие колебания.
Свой вариант задания параметров:
Исходные данные:
Шаг сетки:
Число узлов сетки:
Формирование вектора правой части системы ОДУ и
вектора начальных условий для применения встроенной функции rkfixed:
График решения:
Минимальное значение функции x(t)
достигается в момент времени 0.5.
Максимальное значение функции x(t)
достигается в момент времени 0.1.
Набор параметров подобран таким образом, что
затухающие колебания происходят подобно математическому маятнику -
сопротивление среды останавливает со временем движение груза, происходящее по
гармоническому закону.
Вывод: дифференциальные уравнения второго
порядка - часто используемый способ описания движения. Численное решение этих
дифференциальных уравнений порой единственный способ нахождения закона
движения.
3. Задача
№ 3 (6.2)
3.1 Постановка задачи
Даны две задачи Коши для систем ОДУ 1 порядка с
постоянными коэффициентами на отрезке
[0, 1]
,
,
где A и B - заданные
матрицы, - заданные векторы. Выяснить, какая из задач является жесткой.
ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:
. Составить программу-функцию
нахождения решения системы ОДУ 1 порядка с постоянными коэффициентами по явному
методу Эйлера. Используя составленную программу, решить обе задачи с шагом h=0.01.
Определить, для какой из задач явный метод неустойчив при данном шаге h.
. Используя встроенную функцию eigenvals(M) (M - матрица)
пакета MATHCAD для
нахождения собственных чисел матриц A и B, найти
коэффициенты жесткости обеих систем. Какая из задач является жесткой?
. Для жесткой задачи теоретически
оценить шаг h*, при
котором явный метод Эйлера будет устойчив (см. ПРИЛОЖЕНИЕ C).
. Составить программу-функцию
нахождения решения системы ОДУ 1 порядка с постоянными коэффициентами по
неявному методу Эйлера. Используя составленную программу, найти решение жесткой
задачи с шагом h=0.01. Построить графики компонент
полученного решения.
. Для жесткой задачи
экспериментально подобрать шаг h, при котором графики компонент
решения, полученного по явному методу Эйлера, визуально совпадают с графиками
компонент решения, полученного по неявному методу с шагом h=0.01.
Сравнить найденное значение шага
с шагом h*. Объяснить
различие поведения явного и неявного методов Эйлера при решении жесткой задачи.
УКАЗАНИЕ. В п. 4 для решения системы
линейных уравнений удобно использовать встроенную функцию lsolve пакета MATHCAD.
.2 Исходные данные
|
N
|
A
|
|
B
|
|
|
6.2
|
-17.359 -0.573 5.366 -21.351
|
-64.712 -85.344 -128.964 -170.918
|
1 0
|
3.3 Решение поставленной задачи.
Начальные условия:
Концы отрезка:
Опишем функцию нахождения решения системы ОДУ 1
порядка c постоянными коэффициентами по явному методу Эйлера.
Описанная программа-функция возвращает таблицу
решений, первый столбец которой - это значения аргумента в узлах равномерной
сетки, а остальные столбцы - соответствующие значения компонент приближенного
решения.
Шаг:
Зададим векторы правых частей систем уравнений:
Графики решений для первой системы:
Графики решений для второй системы:
Определим, для какой из задач явный метод
неустойчив при шаге h = 0.01. Найдем собственные числа матриц.
Максимальные и минимальные собственные числа
матриц А и В:
Условие устойчивости выполняется для матрицы А
(usА > h), но не выполняется для матрицы В (usВ < h). Следовательно,
явный метод Эйлера неустойчив для решения системы, описанной матрицей В.
Определим, какая из систем является жесткой:
Число жесткости системы gA
мало (т. е. собственные числа матрицы А незначительно отличаются друг от
друга), поэтому система не жесткая.
Число жесткости системы gB
велико (т. е. собственные числа матрицы В значительно отличаются друг от
друга), поэтому система жесткая.
Определим, при каком шаге явный метод Эйлера
будет устойчив при решении жесткой системы:
задача коши погрешность
дифференциальный
Графики решений для первой и второй компоненты
системы B:
Как видно из графиков решений, явный метод
Эйлера устойчив с шагом hz
= 0.0028 . Условие устойчивости usB>hz
(8.496*10-3 >0.0028) выполняется.
Найдем решение жесткой задачи по неявному методу
Эйлера.
В качестве параметров она принимает матрицу М
системы, вектор начальных условий Vo начало to , конец
отрезка интегрирования T и число узлов равномерной сетки N:
Для оценки результатов решения будем
использовать встроенную функцию для решения жёстких систем stiffr. Для её
применения необходима матрица Якоби:
Графики решений для первой и второй компоненты
системы:
Найдем такое значение шага H для решения жесткой
системы по явному методу Эйлера, что результаты решения будут визуально
совпадать с решением, полученным неявным методом Эйлера с шагом h = 0.01:
Вывод: явный метод Эйлера 1-го порядка точности
дает приближённое решение систем ОДУ с постоянными коэффициентами. При решении
жестких систем ОДУ, метод может быть неустойчив при достаточно большом шаге
вычислений. При уменьшении шага вычислений метод будет устойчив, но это требует
дополнительных (на некотором промежутке лишних) вычислений. Устойчивое решение,
получаемое при решении жёсткой системы уравнений неявным методом, требует в
несколько десятков раз меньше итераций, чем решение, полученное по явному
методу Эйлера.
Заключение
В результате выполнения данной курсовой работы
было реализовано решение задачи Коши с использованием пакета MATHCAD.
При решении различных уравнений были изучены
встроенные функции пакета MATHCAD,
а так же запрограммированы пользовательские функции, позволяющие реализовать
иные методы решения обыкновенных дифференциальных уравнений первого порядка, а
также обыкновенных дифференциальных уравнений второго порядка.
В ходе работы были определены погрешности
решений используемых методов, найдены способы увеличения точности получаемых
результатов.
Так же были построены графики, демонстрирующие
последовательные приближения к искомым решениям.
Таким образом, задание выполнено в полном
объеме.
Список литературы
1)Амосов
А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. М.:
Высшая школа, 1994.
)Сайт
http://www.nsu.ru/matlab/Exponenta_RU/educat/class/courses/vvm/theme_10/theory.asp.htm
)Сайт
http://www.apmath.spbu.ru/ru/education/final/question38.pdf
4)
Ю. Ю. Тарасевич. Численные методы на Mathcad’е. Астрахань, 2000.