Методические указания по МДК.01.04 «Математические методы разработки алгоритмов» (стр. 5 )

Подпись: f(x)

Подпись: f(a) – правая формула (формула правых прямоугольников).

Подпись: f(b) Подпись: f(a) Подпись: f(x) – центральная формула (формула средних прямоугольников).

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

Интервал разобьем на п частей . Тогда – шаг. Удобно работать, когда .

(2.4)

(2.4) – формула прямоугольников.

Погрешность имеет порядок h (~)

Метод трапеций

Аппроксимация в этом методе осуществляется полиномом первой степени. Суть метода ясна из рисунка (рис. 8.). Функция аппроксимируется наклонной линией.

За значение интеграла принимается площадь трапеций:

y(n-1)

y(x1)

y'(x1)

y(n-1)(x1)

y(x2)

y'(x2)

y(n-1)(x2)

…

y(xN)

y'(xN)

y(n-1)(xN)

Например, для дифференциального уравнения первого порядка таблица решения будет представлять собой два столбца – x и y.

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

или , i = 1, …, N

Для определения частного решения необходимо задать дополнительные условия, которые позволят вычислить константы интегрирования. Причем, таких условий должно быть ровно n. Для уравнений первого порядка – одно, для второго – 2 и т. д. В зависимости от способа их задания при решении дифференциальных уравнений существуют два типа задач:

Задача Коши (начальная задача)

Необходимо найти такое частное решение дифференциального уравнения, которое удовлетворяет определенным начальным условиям, заданным в одной точке:

то есть, задано определенное значение независимой переменной (х0) и значение функции и всех ее производных вплоть до порядка (n-1) в этой точке. Эта точка (х0) называется начальной. Например, если решается ДУ 1-го порядка, то начальные условия выражаются в виде пары чисел (x0, y0): .

Краевая задача

В этом случае известны значения функции и (или) ее производных в более чем одной точке, например, в начальный и конечный момент времени, и необходимо найти частное решение дифференциального уравнения между этими точками. Сами дополнительные условия в этом случае называются краевыми (граничными) условиями. Естественно, что краевая задача может решаться для ОДУ не ниже 2-го порядка. Ниже приведен пример ОДУ второго порядка с граничными условиями (заданы значения функции в двух различных точках):

Одношаговые методы решения задачи Коши ОДУ первого порядка

Рассмотрим некоторые одношаговые численные методы решения задачи Коши (начальной задачи) обыкновенных дифференциальных уравнений первого порядка. В одношаговых методах решение в новой точке получается из одного решения в предыдущей точке. Запишем дифференциальное уравнение в общем виде, разрешенном относительно производной (правая часть уравнения не зависит от первой производной):

(3.1)

Необходимо найти значения функции y в заданных точках сетки , если известны начальные значения , где есть значение функции y(x) в начальной точке x0.

Преобразуем уравнение (3.1) умножением на dx

и проинтегрируем левую и правую части между i-ым и i+1-ым узлами сетки:

(3.2)

Мы получили выражение для построения решения в i+1 узле интегрирования через значения x и y в i-ом узле сетки. Сложность, однако, заключается в том, что интеграл в правой части есть интеграл от неявно заданной функции, нахождение которого в аналитическом виде в общем случае невозможно. Численные методы решения ОДУ различным способом аппроксимируют (приближают) значение этого интеграла для построения формул численного интегрирования ОДУ.

Из множества разработанных для решения ОДУ первого порядка одношаговых методов рассмотрим методы Эйлера и Рунге – Кутты. Они достаточно просты и дают начальное представление о подходах к решению данной задачи в рамках численного решения задачи Коши.

Метод Эйлера

Исторически первым и наиболее простым способом численного решения задачи Коши для ОДУ первого порядка является метод Эйлера. В его основе лежит аппроксимация производной отношением конечных приращений зависимой (y) и независимой (x) переменных между узлами равномерной сетки:

где yi+1 это искомое значение функции в точке xi+1.

Если теперь преобразовать это уравнение, и учесть равномерность сетки интегрирования, то получится итерационная формула Эйлера, по которой можно вычислить yi+1 , если известно yi в точке хi:

(3.3)

Сравнивая формулу Эйлера с общим выражением (3.2), полученным ранее, видно, что для приближенного вычисления интеграла в методе Эйлера используется простейшая формула интегрирования – формула левых прямоугольников.

Графическая интерпретация метода Эйлера также не представляет затруднений (рис. 3.1). Действительно, исходя из вида решаемого уравнения (3.1) следует, что значение есть значение производной функции y(x) в точке x = xi , и, таким образом, равно тангенсу угла наклона каcательной, проведенной к графику функции y(x) в точке x = xi.

Рис. 3.1

Из прямоугольного треугольника на рис. 3.1 можно найти:

откуда и получается формула Эйлера. Таким образом, суть метода Эйлера заключается в замене функции y(x) на отрезке интегрирования прямой линией, касательной к графику в точке x=xi. Если искомая функция сильно отличается от линейной на отрезке интегрирования, то погрешность вычисления будет значительной. Ошибка метода Эйлера прямо пропорциональна шагу интегрирования:

Ошибка ~ h

Процесс вычислений строится следующим образом. При заданных начальных условиях x0 и y0 можно вычислить

Таким образом, строится таблица значений функции y(x) с определенным шагом (h) по x на отрезке [x0, xN]. Ошибка в определении значения y(xi) при этом будет тем меньше, чем меньше выбрана длина шага h (что определяется точностью формулы интегрирования).

Метод Эйлера – Коши с итерациями

Точность метода Эйлера можно повысить, если воспользоваться для аппроксимации интеграла более точной формулой интегрирования –формулой трапеций:

(3.4)

Данная формула оказывается неявной относительно yi+1 (это значение есть и в левой и в правой части выражения), то есть является уравнением относительно yi+1, решать которое можно, например, численно, применяя какой-либо итерационный метод (в таком виде его можно рассматривать как итерационную формулу метода простой итерации). Однако можно поступить иначе и приблизительно вычислить значение функции в узле i+1 с помощью обычной формулы Эйлера:

которое затем использовать при вычислении по (3.4).

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

(3.5)

Благодаря более точной формуле интегрирования, погрешность метода пропорциональна уже квадрату шага интегрирования:

Ошибка ~ h2

Методы Рунге – Кутты

Дальнейшее улучшение точности решения ОДУ первого порядка возможно за счет увеличения точности приближенного вычисления интеграла в выражении (3.2).

Мы уже видели, какое преимущество дает переход от интегрирования по формуле прямоугольников (метод Эйлера) к использованию формулы трапеций (метод Эйлера – Коши с итерациями) при аппроксимации этого интеграла.

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

В формуле Симпсона для приближенного вычисления определенного интеграла используются значения подынтегрального выражения в трех точках. В интеграле их всего две, поэтому введем дополнительную точку в середине отрезка [xi+1 xi]

тогда (3.2) можно переписать так:

откуда получаем

Полученное выражение является неявным, так как в правой части содержатся еще не определенные значения функции yi+h/2 и yi+1. Чтобы воспользоваться этой формулой, надо использовать некоторое приближение для вычисления этих значений

При использовании различных методов приближенного вычисления этих величин получаются выражения для методов Рунге-Кутты различного порядка точности.

Формулы метода Рунге-Кутты третьего порядка – РК3 (погрешность порядка h3):

(3.6)

где

Формулы метода Рунге-Кутты четвертого порядка – РК4 (погрешность порядка h4):

, (3.7)

где

Метода Рунге – Кутты третьего и четвертого порядков требуют на каждом шаге трех и четырех вычислений функции соответственно, но являются весьма точными.

Оценка погрешности одношаговых методов

Абсолютная погрешность одношаговых методов d по сравнению с неизвестным точным решением может быть оценена по правилу Рунге

где – решение уравнения с шагом h, – решение уравнения с шагом h/2, p – порядок метода (для метода Эйлера p = 1, для метода Эйлера – Коши с итерациями p = 2, для метода Рунге – Кутты четвертого порядка p = 4).

Порядок выполнения лабораторной работы

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

ПРИМЕРНЫЕ ВОПРОСЫ НА ЗАЩИТЕ РАБОТЫ

1. Что значит – решить задачу Коши для дифференциальных уравнений первого порядка?

2. Графическая интерпретация численного решения дифференциального уравнения.

3. Какие существуют методы решения дифференциального уравнения в зависимости от формы представления решения?

4. В чем заключается суть метода Эйлера?

5. Применение каких формул позволяет получить значения искомой функции по методу Эйлера?

6. Графическая интерпретация метода Эйлера и метода Эйлера – Коши с итерациями. В чем отличие?

7. В чем заключается суть методов Рунге-Кутты?

8. Как определить количество верных цифр в числе, являющемся решением дифференциального уравнения методами Эйлера, Эйлера – Коши с итерациями, Рунге-Кутты?

Варианты заданий

1. Решить методами Эйлера и Рунге – Кутты четвертого порядка на отрезке [0, 0.6] дифференциальное уравнение

при с шагом .

Используя формулу Рунге, оценить погрешность решения.

Построить графики решений.

2. Решить методами Эйлера – Коши с итерациями и Рунге – Кутты третьего порядка дифференциальное уравнение

Построить графики решений.

Полученное численное решение сравнить с точным решением

3. Методами Эйлера – Коши с итерациями и Рунге – Кутты четвертого порядка на отрезке [0, 1] решить с заданной точностью (вводится с клавиатуры) дифференциальное уравнение

при .

Построить графики решений.

4. Решить методами Эйлера и Рунге – Кутты четвертого порядка на отрезке [0, 3] с шагом дифференциальное уравнение

при начальном условии .

Используя формулу Рунге, оценить погрешность решения.

Построить графики решений.

5. Методами Эйлера – Коши с итерациями и Рунге – Кутты четвертого порядка на отрезке [0, 1] решить с заданной точностью (вводится с клавиатуры) дифференциальное уравнение

при начальном условии .

Построить графики решений.

6. Решить методами Эйлера и Рунге – Кутты третьего порядка на отрезке [-1, 1] с шагом дифференциальное уравнение

при начальном условии .

Используя формулу Рунге, оценить погрешность решения.

Построить графики решений.

7. Решить методами Эйлера и Рунге – Кутты третьего порядка дифференциальное уравнение

Построить графики решений.

Полученное численное решение сравнить с точным решением

8. Решить методами Эйлера и Рунге – Кутты четвертого порядка на отрезке [0, 2] с шагом дифференциальное уравнение

при начальном условии .

Используя формулу Рунге, оценить погрешности решений.

Построить графики решений.

9. Решить методами Эйлера – Коши с итерациями и Рунге – Кутты четвертого порядка дифференциальное уравнение

Построить графики решений.

Полученное численное решение сравнить с точным решением

10. Методами Эйлера и Рунге – Кутты четвертого порядка на отрезке [0, 3] найти с заданной точностью (вводится с клавиатуры) решение дифференциального уравнения

при .

Построить графики решений.

11. Решить дифференциальное уравнение

при

на отрезке [0, 1] с шагом :

– методом Рунге – Кутты четвертого порядка;

– методом Эйлера – Коши с итерациями.

Используя формулу Рунге, оценить погрешность решения.

Построить графики решений.

12. Решить методами Эйлера и Рунге – Кутты четвертого порядка дифференциальное уравнение

Построить графики решений.

Полученное численное решение сравнить с точным решением

13. Решить методами Эйлера и Эйлера – Коши с итерациями на отрезке [0, 1] с шагом дифференциальное уравнение

при начальном условии .

Используя формулу Рунге, оценить погрешность решения.

Построить графики решений.

14. Методами Эйлера – Коши с итерациями и Рунге – Кутты третьего порядка на отрезке [0, 3] найти с заданной точностью (вводится с клавиатуры) решение дифференциального уравнения

при .

Построить графики решений.

15. Решить дифференциальное уравнение

при

на отрезке [0, 1] с шагом :

– методом Рунге – Кутты четвертого порядка;

– методом Эйлера.

Используя формулу Рунге, оценить погрешность решения.

Построить графики решений.

16. Решить методами Эйлера – Коши с итерациями и Рунге – Кутты четвертого порядка дифференциальное уравнение

Построить графики решений.

Полученное численное решение сравнить с точным решением

17. Методами Эйлера и Рунге – Кутты четвертого порядка на отрезке [1, 2] решить с заданной точностью (вводится с клавиатуры) дифференциальное уравнение

при начальном условии .

Построить графики решений.

18. Решить методами Рунге – Кутты третьего и четвертого порядков на отрезке [0, 1.5] с шагом дифференциальное уравнение

при начальном условии .

Используя формулу Рунге, оценить погрешность решения.

Построить графики решений.

19. Решить методами Эйлера – Коши с итерациями и Рунге – Кутты третьего порядка дифференциальное уравнение

, .

Построить графики решений.

Полученное численное решение сравнить с точным решением

20. Методами Рунге – Кутты третьего и четвертого порядков на отрезке [0, 1] решить с заданной точностью (вводится с клавиатуры) дифференциальное уравнение

при начальном условии .

Построить графики решений.

Лабораторная работа № 4

Решение систем линейных алгебраических уравнений

Цель работы – изучение прямых и итерационных методов численного решения систем линейных алгебраических уравнений (СЛАУ), приобретение навыков программирования реализации этих методов на языке С.

Методические указания

Рассмотрим систему n линейных алгебраических уравнений с n неизвестными:

В ней aij – коэффициенты при неизвестных xj. Решением этой системы называется такой набор значений неизвестных xj, который удовлетворяет системе.

Коэффициенты aij можно записать в виде матрицы (таблицы):

, правую часть системы в виде вектора , а неизвестные в виде вектора . Тогда систему можно записать в виде матрично-векторного уравнения .

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

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

Метод Гаусса

Этот метод относится к прямым методам решения линейных систем. Он основан на приведении матрицы системы к треугольному виду путем последовательного исключения неизвестных из уравнений системы (прямой ход метода Гаусса) и последующем решении этой треугольной системы, начиная с последнего уравнения (обратный ход метода Гаусса).

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

Этот процесс называется прямым ходом Гаусса и продолжается до тех пор, пока в левой части последнего (n-го) уравнения не останется лишь один член с неизвестным .

Если на каком-то этапе этого процесса оказывается, что очередной исключаемой переменной уже нет ни в одном из последующих уравнений, то матрица системы является вырожденной, и метод Гаусса в этом случае неприменим.

Обратный ход метода Гаусса состоит в последовательном вычислении неизвестных. Решая последнее уравнение, находят единственное неизвестное . Далее, используя это значение, из предыдущего уравнения вычисляют и т. д. Последним находят из первого уравнения.

Рассмотрим применение метода Гаусса для системы из трех уравнений:

(1)

Для исключения из второго уравнения прибавим к нему первое, умноженное на . Затем, умножив первое уравнение на и прибавив результат к третьему уравнению, также исключим из него . Получим равносильную систему уравнений вида:

(2)

Теперь из третьего уравнения системы (2) нужно исключить . Для этого умножим второе уравнение на и прибавим результат к третьему. Получим:

(3)

Матрица системы (3) имеет треугольный вид. На этом завершается прямой ход метода Гаусса.

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

Обратный ход начинается с решения третьего уравнения системы (3):

Используя это значение, можно найти из второго уравнения, а затем из первого:

Аналогично строится вычислительный алгоритм для линейной системы с другим числом неизвестных.

Краткое описание метода Гаусса

Метод Гаусса относится к точным методам решения систем линейных уравнений. Алгоритм этого метода состоит из двух основных этапов: прямого хода и обратного хода.

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

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

, (1)

где – коэффициенты -й строки матрицы после деления; –коэффициенты -й строки исходной матрицы; – номер столбца, ; – правая часть -го уравнения исходной матрицы; – правая часть -го уравнения после преобразования.

Затем производится преобразование всех коэффициентов и правых частей уравнений с по -е с помощью выражений:

(2)

где – номер уравнения ; – номер столбца в уравнении .

После завершения этих преобразований алгоритм повторяется для следующего . После завершения -гo шага начинается обратный ход – расчет неизвестных по выражениям:

(3)

Краткое описание метода LU разложения

Квадратную матрицу можно представить в виде произведения двух треугольных матриц: , где

–нижняя треугольная матрица,

– верхняя треугольная матрица.

С помощью LU разложения матрицы можно получить обратную матрицу и решить систему в виде .

Вычисление коэффициентов треугольных матриц:

(4)

где – соответствующие коэффициенты матрицы .

Коэффициенты обратной матрицы рассчитываются с помощью коэффициентов треугольных матриц и :

. (5)

Итерационные методы

Для применения итерационных методов необходимо предварительно исходную систему уравнений (1) привести к виду:

, , , ,

Этот вид получается, если из первого уравнения выразить x1, из второго x2 и т. д.:

Пусть – некий произвольно задаваемый вектор начального приближения к решению системы. Тогда для нахождения последующих приближений , где m – номер итерации, а , можно применить один из следующих известных методов: метод простой итерации, метод Гаусса-Зейделя, метод верхней релаксации.

1) Метод простой итерации (метод Якоби)

В этом методе коэффициенты вектора рассчитываются по формуле:

2) Метод Гаусса-Зейделя

В этом методе коэффициенты вектора рассчитываются по формуле:

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

При этом хотя бы для одного уравнения неравенство должно выполняться строго.

Эти условия являются достаточными, но не необходимыми, т. е. для некоторых систем итерационный процесс сходится и при нарушении этих условий.

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

Пример 1

Найдем решение следующей линейной системы методом Гаусса:

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

Но мы для простоты понимания проделаем это в два этапа. Сначала сделаем коэффициенты перед переменной x1 во всех уравнениях равными единице, поделив каждое уравнение на коэффициент, стоящий перед этой переменной. Т. е. поделив первое на 2, второе на 2, а третье на 4.

Поместив слева схему производимых действий, запишем полученную систему, эквивалентную исходной:

Теперь избавимся от переменной x1 во втором и третьем уравнениях, вычтя из них первое:

Теперь нам нужно с помощью второго уравнения избавиться от переменной x2 в третьем уравнении. Сделаем это тоже в два этапа. Сначала поделим второе уравнение на 0.5, а третье уравнение на -1.75. Получим систему:

Далее преобразуем третье уравнение, вычтя из него второе:

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

Сделаем проверку, подставив полученные значения неизвестных в левые части уравнений системы, чтобы убедиться в выполнении условий:

Пример 2

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

Имеем систему:

Поставим в ней третье уравнение на первое место, второе уравнение – на третье место, а на втором месте запишем разность второго уравнения и удвоенного первого:

Для такой системы достаточные условия сходимости уже выполняются.

Выразим неизвестные из уравнений, как это предлагается в методе:

Пусть первым приближением решения будет вектор . Тогда следующее приближение рассчитаем по полученным формулам:

Таким образом, имеем следующее приближение решения:

Следующее приближение:

Таким образом, . Далее

Таким образом, .

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