Ирина Эланс

Автор который поможет с любыми образовательными и учебными заданиями

Численное решение задачи Коши

Министерство образования и науки Российской Федерации

ПЕНЗЕНСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ

Кафедра «Системы автоматизированного проектирования»

Пояснительная записка

к курсовой работе по дисциплине

«Вычислительная математика»

на тему: "Численное решение задачи Коши".

Автор работы: Аминева Н.

Специальность «Математическое обеспечение и администрирование информационных систем»

Группа 09ВА1

Руководитель работы Гудков П.А.

Работа защищена «__» ____ 2012 г. Оценка ______________

Пенза 2012 г.

Аннотация.

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

Пояснительная записка содержит 27 страниц, 15 графиков.

Оглавление.

Аннотация. 2

Оглавление. 3

1. Задача № 1 (1.4) 4

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

1.2. Исходные данные. 4

1.3. Решение поставленной задачи. 4

2. Задача № 2 (2.2) 11

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

2.2. Исходные данные. 11

2.3. Решение поставленной задачи. 12

3. Задача № 3 (6.2) 18

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

3.2. Исходные данные. 18

3.3. Решение поставленной задачи. 19

Заключение. 26

Список литературы. 27

1. Задача № 1 (1.4)

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

Найти приближенное решение задачи Коши для обыкновенного дифференциального

уравнения (ОДУ) 1 порядка

(1)

и оценить погрешность решения задачи.

Порядок решения задачи:

1. Задать исходные данные: функцию f правой части, начальное значение .

2. Используя функцию eyler (см. ПРИЛОЖЕНИЕ B), найти приближенное решение задачи Коши с шагом h=0.1 по явному методу Эйлера.

3. Используя встроенную функцию rkfixed пакета MATHCAD, найти приближенное решение задачи Коши с шагом h=0.1 по методу Рунге-Кутты 4 порядка точности (см. ПРИЛОЖЕНИЕ B).

4. Найти решение задачи Коши аналитически.

5. Построить таблицы значений приближенных и точного решений. На одном чертеже построить графики приближенных и точного решений.

6. Оценить погрешность приближенных решений двумя способами:

a) по формуле ; здесьи - значения точного и приближенного решений в узлах сетки , i=1,..N;

b) по правилу Рунге (по правилу двойного пересчета) (см. ПРИЛОЖЕНИЕ C).

7. Выяснить, при каком значении шага h=h* решение, полученное по методуЭйлера, будет иметь такую же погрешность (см. п. 6а), как решение, полученное с помощью метода Рунге-Кутты с шагом h=0.1.

УКАЗАНИЕ. В п. 7 рекомендуется провести серию вычислений решения по методу Эйлера, дробя шаг h пополам.

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

f(t,y)

1.4

Решение поставленной задачи.

Задача Коши: y’(t)=, t0=0, T=1, y0=1.

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

Начальное значение:

Концы отрезка:

Шаг сетки:

Число узлов сетки:

Функция, реализующая явный метод Эйлера, возвращает вектор решения:

Входные параметры:

f - функция правой части;

y0 - начальное значение;

t0 - начальная точка отрезка;

h - шаг сетки;

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

Приближенное решение задачи Коши с шагом h=0.1 по методу Рунге-Кутты 4 порядка точности с помощью встроенной функции rkfixed пакета MATHCAD.

Функция rkfixed возвращает матрицу, первый столбец которой содержит узлы сетки, а второй – приближенное решение в этих узлах.

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

= ,

По методу вариации произвольной постоянной заменим постоянную С на функцию C(t) и решим неоднородное уравнение:

Подставляем в исходное уравнение:

Решение в MathCad:

Решения, полученные различными способами:

Метод Эйлера:

Метод Рунге-Кутты:

Точное решение:

Графики приближенных и точного решений:

Рассчитаем погрешность полученных приближенных решений:

Погрешность метода Эйлера:

Вычисление погрешности по правилу Рунге:

Вычисление приближенных решений с шагом h/2:

Вычисление погрешностей:

Значение погрешностей:

Проведём серию вычислений решения по методу Эйлера, дробя шаг h пополам.

Первая итерация:

Вторая итерация:

Третья итерация:

И т.д.

Девятая итерация:

При значении шага hd=h/1024=0.1/1024=0,0000977 решение, полученное по методу Эйлера, будет иметь примерно такую же погрешность, как и решение, полученное с помощью метода Рунге-Кутты с шагом h=0.1.

Погрешность решения по методу Эйлера с шагом hd=0,0000977 :

Погрешность решения по методу Рунге-Кутты с шагом h=0.1 :

Вывод: явный метод Эйлера - это численный метод 1-го порядка точности. Метод Рунге-Кутты - это метод 4-го порядка точности. Это означает, что при одном и том же значении шага, метод Рунге-Кутты даёт более точное значение. Поэтому погрешности методов сильно (на несколько порядков) отличаются. В рассмотренном выше примере с помощью метода Рунге-Кутты было получено решение, которое совпадает с решением, полученным аналитическим путём.

Задача № 2 (2.2)

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

Задача Коши для ОДУ 2 порядка

описывает движение груза массы m, подвешенного к концу пружины. Здесь x(t) – смещение груза от положения равновесия, H – константа, характеризующая силу сопротивления среды, k –коэффициент упругости пружины, f(t) – внешняя сила. Начальные условия: – смещение груза в начальный момент времени t=0, – скорость груза в начальный момент времени. Промоделировать движение груза на временном отрезке [0,T] при заданных в индивидуальном варианте трех наборах (I, II, III) значений параметров задачи. Для каждого набора по найденной таблице (или графику) решения задачи определить максимальное и минимальное значения

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

ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:

1. Заменить исходную задачу эквивалентной задачей Коши для системы ОДУ 1 порядка:

(2)

2. Для каждого варианта выбора параметров решить задачу (2) с помощью метода Рунге-Кутты 4 порядка точности с шагом h=0.1.

3. Для каждого варианта выбора параметров построить график найденного решения. Сравнить характер движения груза и дать интерпретацию полученного движения.

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

УКАЗАНИЕ. В п. 2 использовать встроенную функцию rkfixed пакета MATHCAD (см. ПРИЛОЖЕНИЕ B).

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

f(t)

2.2

III

0.5

tsin(t)

-10

-50

Решение поставленной задачи.

1 набор.

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

Шаг сетки:

Число узлов сетки:

Формирование вектора правой части системы ОДУ и вектора начальных условий для применения встроенной функции rkfixed:

График решения:

Найдем максимальное и минимальное значения функции x(t):

Минимальное значение достигается в момент времени 18.4.

Максимальное значение достигается в момент времени 15.3.

Из графика видно, что при данном наборе значений груз совершает незатухающие колебания.

2 набор.

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

Шаг сетки:

Число узлов сетки:

График решения:

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

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

При данном наборе значений происходит затухание колебаний – груз останавливается.

3 набор.

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

Шаг сетки:

Число узлов сетки:

График решения:

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

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

При данном наборе значений, как и в первом случае, груз совершает незатухающие колебания.

Свой вариант задания параметров:

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

Шаг сетки:

Число узлов сетки:

График решения:

Минимальное значение функции x(t) достигается в момент времени 0.5.

Максимальное значение функции x(t) достигается в момент времени 0.1.

Набор параметров подобран таким образом, что затухающие колебания происходят подобно математическому маятнику - сопротивление среды останавливает со временем движение груза, происходящее по гармоническому закону.

Вывод: дифференциальные уравнения второго порядка - часто используемый способ описания движения. Численное решение этих дифференциальных уравнений порой единственный способ нахождения закона движения.

Задача № 3 (6.2)

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

Даны две задачи Коши для систем ОДУ 1 порядка с постоянными коэффициентами на отрезке

[0, 1]

где A и B – заданные матрицы, - заданные векторы. Выяснить, какая из задач является жесткой.

ПОРЯДОК РЕШЕНИЯ ЗАДАЧИ:

1. Составить программу-функцию нахождения решения системы ОДУ 1 порядка с постоянными коэффициентами по явному методу Эйлера. Используя составленную программу, решить обе задачи с шагом h=0.01. Определить, для какой из задач явный метод неустойчив при данном шаге h.

2. Используя встроенную функцию eigenvals(M) (M – матрица) пакета MATHCAD для нахождения собственных чисел матриц A и B, найти коэффициенты жесткости обеих систем. Какая из задач является жесткой?

3. Для жесткой задачи теоретически оценить шаг h*, при котором явный метод Эйлера будет устойчив (см. ПРИЛОЖЕНИЕ C).

4. Составить программу-функцию нахождения решения системы ОДУ 1 порядка с постоянными коэффициентами по неявному методу Эйлера. Используя составленную программу, найти решение жесткой задачи с шагом h=0.01. Построить графики компонент полученного решения.

5. Для жесткой задачи экспериментально подобрать шаг h, при котором графики компонент решения, полученного по явному методу Эйлера, визуально совпадают с графиками компонент решения, полученного по неявному методу с шагом h=0.01. Сравнить найденное значение шага

с шагом h*. Объяснить различие поведения явного и неявного методов Эйлера при решении жесткой задачи.

УКАЗАНИЕ. В п. 4 для решения системы линейных уравнений удобно использовать встроенную функцию lsolve пакета MATHCAD.

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

6.2

-17.359 -0.573

5.366 -21.351

₂

₁

-64.712 -85.344

-128.964 -170.918

₁

₀

Решение поставленной задачи.

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

Концы отрезка:

Опишем функцию нахождения решения системы ОДУ 1 порядка c постоянными коэффициентами по явному методу Эйлера.

Описанная программа-функция возвращает таблицу решений, первый столбец которой - это значения аргумента в узлах равномерной сетки, а остальные столбцы - соответствующие значения компонент приближенного решения.
Шаг:

Зададим векторы правых частей систем уравнений:

Графики решений для первой системы:

Графики решений для второй системы:

Определим, для какой из задач явный метод неустойчив при шаге h = 0.01. Найдем собственные числа матриц.

Максимальные и минимальные собственные числа матриц А и В:

Условие устойчивости выполняется для матрицы А (usА > h), но не выполняется для матрицы В (usВ < h). Следовательно, явный метод Эйлера неустойчив для решения системы, описанной матрицей В.

Определим, какая из систем является жесткой:

Число жесткости системы gA мало (т. е. собственные числа матрицы А незначительно отличаются друг от друга), поэтому система не жесткая.

Число жесткости системы gB велико (т. е. собственные числа матрицы В значительно отличаются друг от друга), поэтому система жесткая.

Определим, при каком шаге явный метод Эйлера будет устойчив при решении жесткой системы:

Графики решений для первой и второй компоненты системы B:

Как видно из графиков решений, явный метод Эйлера устойчив с шагом hz = 0.0028 . Условие устойчивости usB>hz (8.496*10^-3 >0.0028) выполняется.

Найдем решение жесткой задачи по неявному методу Эйлера.

Опишем функцию нахождения решения системы ОДУ 1 порядка c постоянными коэффициентами по неявному методу Эйлера.

В качестве параметров она принимает матрицу М системы, вектор начальных условий V_o начало t_o , конец отрезка интегрирования T и число узлов равномерной сетки N:

Для оценки результатов решения будем использовать встроенную функцию для решения жёстких систем stiffr. Для её применения необходима матрица Якоби:

Графики решений для первой и второй компоненты системы: