• Регистрация
Редактор-сообщества-Экспонента
Редактор-сообщества-Экспонента+26.70
н/д
  • Написать
  • Подписаться

Курс "Введение в вычислительную математику"

Математика и статистика 
21.12.2019

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

1. Теория погрешностей и машинная арифметика.

2. Решение нелинейных уравнений. Методы бисекций, простой итерации, Ньютона.

3. Решение нелинейных уравнений. Обусловленность задачи нахождения корня. Интервал неопределенности.

4. Решение систем линейных алгебраических уравнений. Нормы векторов и матриц. LU-разложение матрицы.

5. Решение систем линейных алгебраических уравнений прямыми методами.

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

7. Приближение функций. Метод наименьших квадратов.

8. Приближение функций. Интерполяция.

9. Приближение функций. Сплайны.

10. Решение задачи Коши одношаговыми методами. 

 

 

Наверх

1. Теория погрешностей и машинная арифметика.

Пусть image2 (217 bytes) - точное значение, 

       image3 (103 bytes) - приближенное значение некоторой величины. 

Абсолютной погрешностью приближенного значения image3 (103 bytes) называется величина image4 (371 bytes) . 

Относительной погрешностью значения image3 (103 bytes) (при image5 (259 bytes) 0) называется величина   image6 (293 bytes)

Так как, значение image2 (217 bytes) как правило неизвестно, чаще получают оценки погрешностей вида: image7 (529 bytes)     image8 (621 bytes)

Величины image9 (286 bytes) и image10 (306 bytes) называют верхними границами (или просто границами) абсолютной и относительной погрешностей. 

 

Значащими цифрами числа  image3 (103 bytes)   называют все цифры в его записи, начиная с первой ненулевой слева.

 

Значащую цифру числа image3 (103 bytes)  называют верной, если абсолютная погрешность числа не превосходит единицы разряда, соответствующего этой цифре. 

 

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

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

image11 (616 bytes)

Если а и b - ненулевые числа одного знака, то справедливы неравенства

image12 (403 bytes),      image13 (410 bytes) , 

где image14 (895 bytes) ,     image15 (404 bytes)

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

если image16 (295 bytes) и image17 (298 bytes) , то   image18 (585 bytes), image19 (626 bytes).

 

Пусть image20 (494 bytes) - дифференцируемая в области G функция переменных, вычисление которой производится при приближенно заданных значениях аргументов   image21 (255 bytes) . Тогда для абсолютной погрешности функции image22 (301 bytes) справедлива следующая оценка

image23 (733 bytes) image24 (435 bytes).

Здесь [x, x*] v отрезок, соединяющий точки x и x* =( image25 (255 bytes)

Для относительной погрешности функции справедливо следующее приближенное равенство

image26 (803 bytes) , где image27 (794 bytes)

 

 

Наверх

2. Решение нелинейных уравнений. Методы бисекций, простой итерации, Ньютона.

Пусть рассматривается уравнение image1.gif(983 bytes). Корнем уравнения называется значение image2.gif(860 bytes), при котором image3.gif(994 bytes). Корень image2.gif(860 bytes) называется простым, если 

image4.gif(1004 bytes), в противном случае корень называется кратным. Целое число m называется кратностью корня image2.gif(860 bytes), если image5.gif(1034 bytes) для k=1,2,3-,m-1 и 

image6.gif(1044 bytes).

Постановка задачи вычисления приближенного значения корня с точностью image7.gif(848 bytes): найти такое значения image8.gif(850 bytes), что image9.gif(1022 bytes).

Решение задачи разбивается на два этапа: на первом этапе осуществляют локализацию корней, на втором этапе производят итерационное уточнение корней. На этапе локализации корней находят достаточно узкие отрезки ( или отрезок, если корень единственный), которые содержат один и только один корень уравнения image1.gif(983 bytes). На втором этапе вычисляют приближенное значение корня с заданной точностью. Часто вместо отрезка локализации достаточно указать начальное приближение к корню.

 

Метод бисекции. Пусть [a,b] v отрезок локализации. Предположим, что функция f(x) непрерывна на [a,b] и на концах принимает значения разных знаков image10.gif(1087 bytes)

 

Алгоритм метода бисекции состоит в построении последовательности вложенных отрезков, на концах которых функция принимает значения разных знаков. Каждый последующий отрезок получают делением пополам предыдущего. Опишем один шаг итераций метода. Пусть на k-ом шаге найден отрезок image11.gif(1043 bytes) такой, что image12.gif(1260 bytes). Найдем середину отрезка image13.gif(1160 bytes). Если image14.gif(1031 bytes), то image15.gif(904 bytes)- корень и задача решена. Если нет, то из двух половин отрезка выбираем ту, на концах которой функция имеет противоположные знаки:

image16.gif(1013 bytes), image17.gif(1012 bytes), если image18.gif(1249 bytes)

 

image19.gif(1012 bytes), image20.gif(1017 bytes), если image21.gif(1256 bytes)

 

Критерий окончания итерационного процесса: если длина отрезка локализации меньше 2image7.gif(848 bytes), то итерации прекращают и в качестве значения корня с заданной точностью принимают середину отрезка.

 

Теорема о сходимости метода бисекций. Пусть функция f(x) непрерывна на [a,b] и на концах принимает значения разных знаков image10.gif(1087 bytes).Тогда метод сходится и справедлива оценка погрешности : image22.gif(1492 bytes)

 

Метод Ньютона (метод касательных) . Расчетная формула метода Ньютона имеет вид:

image23.gif(1396 bytes). Геометрически метод Ньютона означает, что следующее приближение к корню image24.gif(929 bytes) есть точка пересечения с осью ОХ

касательной, проведенной к графику функции y=f(x) в точке image25.gif(1142 bytes)

 

Теорема о сходимости метода Ньютона. Пусть image2.gif(860 bytes) - простой корень уравнения image1.gif(983 bytes), в некоторой окрестности которого функция дважды непрерывно дифференцируема. Тогда найдется такая малая image26.gif(855 bytes)- окрестность корня image2.gif(860 bytes), что при произвольном выборе начального приближения image27.gif(905 bytes) из этой окрестности итерационная последовательность метода Ньютона не выходит за пределы окрестности и справедлива оценка 

image28.gif(1281 bytes), где image29.gif(918 bytes), image30.gif(935 bytes).

 

Критерий окончания итерационного процесса. При заданной точности image7.gif(848 bytes)>0 вычисления следует вести до тех пор пока не окажется выполненным неравенство image31.gif(1124 bytes)

 

Как указано в теореме, метод Ньютона обладает локальной сходимостью, то есть областью его сходимости является малая окрестность корня image2.gif(860 bytes). Неудачный выбор может дать расходящуюся итерационную последовательность.

 

Метод простой итерации (метод последовательных повторений). Для применения метода простой итерации следует исходное уравнение image1.gif(983 bytes) преобразовать к виду, удобному для итерации image32.gif(985 bytes). Это преобразование можно выполнить различными способами. Функция image33.gif(943 bytes) называется итерационной функцией. Расчетная формула метода простой итерации имеет вид: image34.gif(1137 bytes)

 

Теорема о сходимости метода простой итерации. Пусть в некоторой image26.gif(855 bytes)- окрестности корня image2.gif(860 bytes)функция image33.gif(943 bytes) дифференцируема и удовлетворяет неравенству image35.gif(1049 bytes), где image36.gif(968 bytes) - постоянная . Тогда независимо от выбора начального приближения из указанной image26.gif(855 bytes)- окрестности итерационная последовательность не выходит из этой окрестности, метод сходится

со скоростью геометрической последовательности и справедлива оценка погрешности:image37.gif(1375 bytes), image38.gif(908 bytes).

 

Критерий окончания итерационного процесса. При заданной точности image7.gif(848 bytes)>0 вычисления следует вести до тех пор пока не окажется выполненным неравенство image39.gif(1246 bytes). Если величина image40.gif(1004 bytes), то можно использовать более простой критерий окончания итераций: image31.gif(1124 bytes).

Ключевой момент в применении метода простой итерации состоит в эквивалентном преобразовании уравнения. Способ, при котором выполнено условие сходимости метода простой итерации, состоит в следующем: исходное уравнение приводится к виду image41.gif(1034 bytes). Предположим дополнительно, что производная image42.gif(870 bytes) знакопостоянна и image43.gif(1064 bytes) на отрезке [a,b]. Тогда при выборе итерационного параметра image44.gif(1159 bytes) метод сходится и значение

image45.gif(1232 bytes) .

 

 

 

Наверх

3. Решение нелинейных уравнений. Обусловленность задачи нахождения корня. Интервал неопределенности.

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

Пусть установлено неравенство image002.gif(382 bytes), где image004.gif(251 bytes) - относительная погрешность входных данных, а image006.gif(256 bytes)- относительная погрешность решения. Тогда image008.gif(195 bytes)- называется абсолютным числом обусловленности задачи. Если же установлено неравенство image010 (382 bytes) между относительными погрешностями данных и решения, то image012 (195 bytes) называют относительным числом обусловленности задачи. 

 

Обычно под числом обусловленности image014.gif (179 bytes) понимают относительное число обусловленности. Если image016.gif(220 bytes), то задачу называют плохо обусловленной.

Обусловленность задачи нахождения корня. Пусть image018.gif(181 bytes) v корень, подлежащий определению. Будем считать, что входными данными для задачи вычисления корня являются значения функции image020.gif(216 bytes). Так как image021gif (216 bytes)v вычисляется приближенно, то обозначим функцию, полученную в действительности через image023.gif(230 bytes). Предположим, что в малой окрестности корня выполняется неравенство: image025.gif (525 bytes). Для близких к image026.gif (181 bytes) значений image028.gif (178 bytes) справедливо равенство image30.gif (497 bytes), следовательно, image032.gif (595 bytes). Это означает, что число обусловленности задачи нахождения корня равно image034.gif (403 bytes). Из последней формулы следует, что чем меньше значение производной функции в точке корня, тем задача хуже обусловлена. В частности, задача нахождения кратного корня имеет число обусловленности - бесконечность. 

 

Интервал неопределенности корня. Если функция image035.gif (216 bytes) непрерывна, то найдется такая малая окрестность image037.gif (293 bytes)  корня image038.gif (181 bytes), имеющая радиус image040.gif (174 bytes), в которой выполнено неравенство image042.gif (291 byes). Это означает, что image044.gif (187 bytes) image045.gif (293 bytes) знак вычисленного значения image046.gif (230 bytes) , вообще говоря не обязан совпадать со знаком image047.gif (216 bytes) и, следовательно, становится невозможным определить, какое именно значение image049.gif (174 bytes) из интервала image050.gif (293 bytes) обращает функцию image052.gif в нуль. Этот интервал называется интервалом неопределенности корня. Очевидно, что радиус интервала неопределенности для простого корня равен image054.gif (631.gif). Аналогично можно показать, что для кратного корня image056.gif (525 bytes). Это означает, что для простого корня радиус интервала неопределенности прямо пропорционален погрешности вычисления функции image058.gif (256 bytes), а для кратного корня image060.gif (357 bytes)

 

Пусть image062.gif (319 bytes). Корень уравнения простой и равен image063.gif (181 bytes) = -0.34729635533861. Тогда image065.gif (302 bytes) и image067.gif (301 bytes). Если image069.gif (340 bytes), то image071.gif (258 bytes). Это означает , что найти корень с точностью меньшей, чем радиус интервала неопределенности, не удастся. 

  

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

image073.gif (694 bytes), где image075.gif (184 bytes) - кратность корня. 

 

Как правило, значение image076.gif (184 bytes) v неизвестно. Используя метод Ньютона, можно узнать кратность корня. Для этого будем задавать значения image077.gif (184 bytes)= 1,2,3 и вычислять значение корня с заданной точностью , одновременно подсчитывая количество итераций для каждого значения image078.gif (184 bytes). При некотором значении image079.gif (184 bytes)число итераций будет минимальным. Это значение image080.gif (184 bytes) и есть кратность корня.

 

 

Наверх

4. Решение систем линейных алгебраических уравнений. Нормы векторов и матриц. LU-разложение матрицы.

Нормы векторов и матриц. Обозначим через image002 (366 bytes)- точное решение системы, а через image004 (420 bytes)- приближенное решение системы. Для количественной характеристики вектора погрешности image006 (246 bytes) введем понятие нормы.

 

Нормой вектора image008 (179 bytes) называется число image010 (225 bytes), удовлетворяющее трем аксиомам: 

1) image012 (294 bytes) причем image010 (225 bytes) = 0 тогда и только тогда, когда image008 (179 bytes)= 0;

2) image016 (392 bytes) для любого вектора image008 (179 bytes) и любого числа image019 (182 bytes);

3) image021 (462 bytes) для любых векторов image008 (179 bytes) и image024 (185 bytes).

Наиболее употребительными являются следующие три нормы:

image026 (447 bytes)  ,  image028 (635 bytes)  ,  image030 (496 bytes) .

Абсолютная и относительная погрешности вектора вводятся с помощью формул:

image032 (404 bytes)  и  image034 (534 bytes).

Нормой матрицы image036 (186 bytes) называется величина image038 (632 bytes). Введенная норма обладает свойствами, аналогичными свойствам нормы вектора:

1) image040 (299 bytes) причем image042 (232 bytes) = 0 тогда и только тогда, когда  A = 0;

2) image044 (396 bytes) для любой матрицы A и любого числа image019 (182 bytes);

3) image047 (467 bytes) для любых матриц A и B;

4) image049 (369 bytes).

Каждой из векторных норм соответствует своя подчиненная норма матрицы:

image051 (602 bytes)  ,  image053 (496 bytes)  ,  image055 (606 bytes) .

В оценках вместо нормы image057 (233 bytes) используется евклидова норма матрицы

image059 (635 bytes)  , так как  image061 (355 bytes) .

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

image063 (414 bytes)  ,  image065 (559 bytes) .

 

Пусть рассматривается система линейных алгебраических уравнений

image067 (1592 bytes) image069 (115 bytes)

В матричной форме записи она имеет вид image071 (233 bytes). Будем предполагать, что матрица системы image036 (186 bytes)задана и является невырожденной. Известно, что в этом случае решение системы существует, единственно и устойчиво по входным данным.

 

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

 

Теорема об оценке погрешности решения по погрешностям входных данных.

Пусть решение системы image073 (233 bytes), а x* - решение системы A*x*=b*, тогда image075 (491 bytes), где image077 (446 bytes) - относительное число обусловленности системы.

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

 

Метод Гаусса. Рассмотрим метод Гаусса (схему единственного деления) решения системы уравнений. Прямой ход состоит из m-1 шагов исключения.

1 Шаг. Исключим неизвестное image079 (182 bytes) из уравнений с номерами i = 2,3,..m. Предположим, что image081 (232 bytes). Будем называть его ведущим элементом 1-го шага.

Найдем величины image083 (292 bytes), i=2,3,...m , называемые множителями 1-го шага. Вычтем последовательно из второго, третьего, ...m vго уравнений системы первое уравнение, умноженное соответственно на image085 (275 bytes). В результате 1-го шага получим эквивалентную систему уравнений:

image087 (1777 bytes)

Аналогично проводятся остальные шаги. Опишем очередной k-ый шаг. Предположим, что ведущий элемент image089 (233 bytes). Вычислим множители к-го шага:image091 (376 bytes), i=k+1,...m и вычтем последовательно из (k+1)-го, ...m v го уравнений системы k-ое уравнение, умноженное соответственно на 

image093 (317 bytes).После (m-1)-го шага исключения получим систему уравнений

image095 (1688 bytes),

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

Обратный ход. Из последнего уравнения системы находим image097 (187 bytes). Подставляя найденное значение image097 (187 bytes) в предпоследнее уравнение, получим image100 (199 bytes). Далее последовательно находим неизвестные image102 (253 bytes).

 

LU  разложение матрицы. Представим матрицу A в виде произведения нижней треугольной матрицы L и верхней треугольной U.

Введем в рассмотрение матрицы 

image104 (991 bytes) и image106 (195 bytes)image108 (187 bytes)

Можно показать, что A = LU. Это и есть разложение матрицы на множители.

 

 

Наверх

5. Решение систем линейных алгебраических уравнений прямыми методами.

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

 

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

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

Метод Гаусса с выбором главного элемента по столбцу (схема частичного выбора). На k-ом шаге прямого хода в качестве ведущего элемента выбирают максимальный по модулю коэффициент image002.gif (205 bytes) при неизвестной image004.gif (186 bytes) в уравнениях с номерами i = k+1, ... , m.Затем уравнение, соответствующее выбранному коэффициенту с номером image006.gif (184 bytes), меняют местами с к-ым уравнением системы для того, чтобы главный элемент занял место коэффициента image008.gif (233 bytes). После этой перестановки исключение проводят как в схеме единственного деления. В этом случае все масштабирующие множители по модулю меньше единицы и схема обладает вычислительной устойчивостью.

 

Метод Холецкого. Если матрица системы является симметричной и положительно определенной, то для решения системы применяют метод Холецкого (метод квадратных корней). В основе метода лежит алгоритм специального LU-разложения матрицы A, в результате чего она приводится к виду A=image010.gif (204 bytes). Если разложение получено, то как и в методе LU-разложения, решение системы сводится к последовательному решению двух систем с треугольными матрицами: image012.gif (226 bytes) и  image014.gif (239 bytes). Для нахождения коэффициентов матрицы L неизвестные коэффициенты матрицы image010.gif (204 bytes) приравнивают соответствующим элементам матрицы A. Затем последовательно находят требуемые коэффициенты по формулам:

image017.gif (277 bytes) , image019.gif (268 bytes) i = 2, 3, ..., m,

image021.gif (332 bytes) , image023.gif (348 bytes) i = 3, 4, ..., m,

...............

image025.gif (476 bytes)  image027.gif (505 bytes)   i = k+1, ... , m.

image029.gif (588 bytes)

 

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

image031.gif (1935 bytes)

Преобразуем первое уравнение системы к виду image033.gif (291 bytes), где image035.gif (263 bytes), image037.gif (269 bytes)

Подставим полученное выражение во второе уравнение системы и преобразуем его к виду image039.gif (297 bytes) и т.д. На i-ом шаге уравнение преобразуется к виду image041.gif (295 bytes), где image043.gif (343 bytes), image045.gif (435 bytes). На m-ом шаге подстановка в последнее уравнение выражения image047.gif (333 bytes) дает возможность определить значение image049.gif (187 bytes):

image051.gif (511 bytes). Значения остальных неизвестных находятся по формулам: image052.gif (295 bytes), i = m-1, m-2, ..., 1.

 

 

Наверх

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

Рассматривается система Ax = b. 

Для применения итерационных методов система должна быть приведена к эквивалентному виду x=Bx+d. Затем выбирается начальное приближение к решению системы уравнений 

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

 

Метод Якоби.

Самый простой способ приведения системы к виду удобному для итерации состоит в следующем: из первого уравнения системы выразим неизвестное  x1, из второго уравнения системы выразим  x2, и т. д. В результате получим систему уравнений с матрицей B, в которой на главной диагонали стоят нулевые элементы, а остальные элементы вычисляются по формулам:

 image010.gif (1014 bytes),i, j = 1, 2, ... n. 

Компоненты вектора d вычисляются по формулам:

 image012.gif (974 bytes),   i = 1, 2, ... n. 

Расчетная формула метода простой итерации имеет вид 

image014.gif (1048 bytes)

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

image016.gif (1380 bytes),  i = 1, 2, ... m.

Критерий окончания итераций в методе Якоби имеет вид: 

image020.gif (1154 bytes) ,   где   image022.gif (1187 bytes)

Если  image024.gif (1008 bytes), то можно применять более простой критерий 

image026.gif (1146 bytes) окончания итераций 

 

Метод Зейделя.

Метод можно рассматривать как модификацию метода Якоби. Основная идея состоит в том, что при вычислении очередного (n+1)-го приближения к неизвестному  xi при i >1 используют уже найденные (n+1)-е приближения к неизвестным  x1, x2, ..., xi - 1, а не n-ое приближение, как в методе Якоби. Расчетная формула метода в покоординатной форме записи выглядит так:

image034.gif (1729 bytes)

 i = 1, 2, ... m.. Условия сходимости и критерий окончания итераций можно взять такими же как в методе Якоби.

 

Пусть матрица системы уравнений A - симметричная и положительно определенная. Тогда при любом выборе начального приближения метод Зейделя сходится. Дополнительных условий на малость нормы некоторой матрицы здесь не накладывается. 

 

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

Если A - симметричная и положительно определенная матрица, то систему уравнений часто приводят к эквивалентному виду: 

x = x - image039.gif (834 bytes)(Ax - b),  image039.gif (834 bytes) - итерационный параметр. 

Расчетная формула метода простой итерации в этом случае имеет вид:

x (n+1) = x n -  image039.gif (834 bytes)(Ax n - b). 

Здесь B = E -  image039.gif (834 bytes)A  и параметр  image039.gif (834 bytes) > 0 выбирают так, чтобы по возможности сделать минимальной величину  image046.gif (929 bytes)

Пусть image048.gif (902 bytes) и image050.gif (902 bytes) - минимальное и максимальное собственные значения матрицы A. Оптимальным является выбор параметра  image052.gif (1099 bytes). В этом случае  image046.gif (929 bytes) принимает минимальное значение равное image055.gif (1216 bytes).

 

 

Наверх

7. Приближение функций. Метод наименьших квадратов.

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

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

 

Постановка задачи приближения функции по методу наименьших квадратов. Пусть функция y=f(x) задана таблицей своих значений: image4.gif(967 bytes), i=0,1,-n. Требуется найти многочлен фиксированной степени m, для которого среднеквадратичное отклонение (СКО) image5.gif(1333 bytes)минимально.

Так как многочлен image6.gif(1201 bytes)определяется своими коэффициентами, то фактически нужно подобрать набор кофициентов image7.gif(964 bytes), минимизирующий функцию image8.gif(1755 bytes).

Используя необходимое условие экстремума, image9.gif(1008 bytes), k=0,1,-m получаем так называемую нормальную систему метода наименьших квадратов: image10.gif(1392 bytes), k=0,1,-m.

Полученная система есть система алгебраических уравнений относительно неизвестных image7.gif(964 bytes). Можно показать, что определитель этой системы отличен от нуля, то есть решение существует и единственно. Однако при высоких степенях m система является плохо обусловленной. Поэтому метод наименьших квадратов применяют для нахождения многочленов, степень которых не выше 5. Решение нормальной системы можно найти, например, методом Гаусса.

Запишем нормальную систему наименьших квадратов для двух простых случаев: m=0 и m=2. При m=0 многочлен примет вид: image11.gif(979 bytes). Для нахождения неизвестного коэффициента image12.gif(872 bytes) имеем уравнение:image13.gif(1194 bytes). Получаем, что коэффициент image14.gif(881 bytes) есть среднее арифметическое значений функции в заданных точках. 

Если же используется многочлен второй степени image15.gif(1105 bytes), то нормальная система уравнений примет вид:

image16.gif(3981 bytes)

Предположим, что функцию f можно с высокой точностью аппроксимировать многочленом image17.gif(929 bytes)некоторой степени m. Если эта степень заранее неизвестна, то возникает проблема выбора оптимальной степени аппроксимирующего многочлена в условиях, когда исходные данные image18.gif(865 bytes) содержат случайные ошибки. Для решения этой задачи можно принять следующий алгоритм: для каждого m=0,1,2,.. вычисляется величина 

image19.gif(1365 bytes). За оптимальное значение степени многочлена следует принять то значение m, начиная с которого величина image20.gif(879 bytes) стабилизируется или начинает возрастать.

 

Определение параметров эмпирической зависимости. Часто из физических соображений следует, что зависимость image21.gif(946 bytes) между величинами хорошо описывается моделью вида image22.gif(1079 bytes), где вид зависимости g известен. Тогда применение критерия наименьших квадратов приводит к задаче определения искомых параметров image7.gif(964 bytes)из условия минимума функции: image23.gif(1436 bytes).

 

Если зависимость от параметров image24.gif(969 bytes) нелинейна, то экстремум функции image23.gif(1436 bytes)ищут методами минимизации функций нескольких переменных.

 

 

Наверх

8. Приближение функций. Интерполяция.

Постановка задачи интерполяции функций.

Пусть функция y = f(x) задана таблицей своих значений: 

image002.gif(1118 bytes), i=0,1,...n. Требуется найти многочлен степени n, такой, что значения функции и многочлена в точках таблицы совпадают:

image004.gif(328 bytes), i=0,1,... n. 

Справедлива теорема о существовании и единственности интерполяционного многочлена.

 

Одна из форм записи интерполяционного многочлена - многочлен Лагранжа:

image006.gif(437 bytes), где 

image008.gif(1296 bytes)

Многочлен image010.gif(229 bytes)представляет собой многочлен степени n , удовлетворяющий условию

image012.gif(447 bytes) .

Таким образом, степень многочлена image014.gif(189 bytes) равна n и при image016.gif(206 bytes) в сумме обращаются в нуль все слагаемые, кроме слагаемого с номером image018.gif(200 bytes), равного image020.gif(181 bytes). Поэтому многочлен Лагранжа является интерполяционным.

 

 

Другая форма записи интерполяционного многочлена - интерполяционный многочлен Ньютона с разделенными разностями. Пусть функция f задана с произвольным шагом и

точки таблицы занумерованы в произвольном порядке. Величины 

image022.gif(534 bytes) называют разделенными разностями первого порядка. Разделенные разности второго порядка определяются формулой:

image024.gif(663 bytes).

Определение разделенной разности порядка image026.gif(210 bytes)таково: 

image028.gif(720 bytes)

Используя разделенные разности, интерполяционный многочлен Ньютона можно записать в следующем виде:

0301.gif (1360 bytes)

0302.gif (1180 bytes)

 

Величину image032.gif(397 bytes)называют погрешностью интерполяции или остаточным членом интерполяции.

 

Оценка погрешности интерполяции.

Если функция n+1 раз на отрезке [a,b] , содержащем узлы интерполяции image034.gif(181 bytes), i=0,1,...n, то для погрешности интерполяции справедлива оценка: 

image036.gif(722 bytes). Здесь image038.gif(520 bytes), image040.gif(372 bytes).

Эта оценка показывает, что для достаточно гладкой функции при фиксированной степени интерполяционного многочлена погрешность интерполяции стремится к нулю не медленнее, чем величина, пропорциональная image042.gif(225 bytes). Этот факт формулируют так: интерполяционный многочлен степени nаппроксимирует функцию с (n+1) порядком точности относительно image044.gif(211 bytes).

 

В практическом плане формула Ньютона обладает преимуществами перед формулой Лагранжа. Предположим, что в необходимо увеличить степень многочлена на единицу, добавив в таблицу еще один узел image046.gif(197 bytes). При использовании формулы Лагранжа это приводит к необходимости вычислять каждое слагаемое заново. Для вычисления image048.gif(239 bytes)достаточно добавить к image050.gif(225 bytes)лишь очередное слагаемое: image052.gif(496 bytes). Если функция f достаточно гладкая, то справедливо приближенное равенство image054.gif(433 bytes). Это равенство можно использовать для практической оценки погрешности интерполяции :image056.gif(376 bytes).

 

 

Наверх

9. Приближение функций. Сплайны.

Глобальная и кусочно-полиномиальная интерполяция. Пусть функция  f(x) интерполируется на отрезке [a,b]. Метод решения этой задачи с помощью единого многочлена image002.gif (225 bytes) для всего отрезка называют глобальной полиномиальной интерполяцией. В вычислительной практике такой поход применяется редко в силу различных причин. Одна из причин v необходимо задать стратегию выбора узлов при интерполяции функции  f многочленами все возрастающей степени n.

Теорема Фабера. Какова бы ни была стратегия выбора узлов интерполяции, найдется непрерывная на отрезке [a,b] функция  f, для которой image004.gif (512 bytes) при image006.gif (215 bytes). Таким образом, теорема Фабера отрицает существование единой для всех непрерывных функций стратегии выбора узлов. Проиллюстрируем сказанное примером.

Предположим, что выбираем равноотстоящие узлы, то есть imge008.gif (252 bytes), i = 0,1,...n, где image010.gif (283 bytes). Покажем, что для функции Рунге такая стратегия является неудачной.

 

На практике чаще используют кусочно-полиномиальную интерполяцию: исходный отрезок разбивается на части и на каждом отрезке малой длины исходная функция заменяется многочленом невысокой степени. Система Mathcad предоставляет возможность аппроксимации двумя важными классами функций: кусочно-линейной и сплайнами.

При кусочно-линейной интерполяции узловые точки соединяются отрезками прямых, то есть через каждые две точки imge012.gif (265 bytes) , imge014.gif (292 bytes) проводится полином первой степени.

Как видно из приведенного примера этот способ приближения также имеет недостаток: в точках "стыка" двух соседних многочленов производная, как правило, имеет разрыв.

Если исходная функция была гладкой и требуется, чтобы и аппроксимирующая функция была гладкой, то кусочно-полиномиальная интерполяция неприемлема. В этом случае применяют сплайны v специальным образом построенные гладкие кусочно-многочленные функции.

Интерполяция сплайнами. Пусть отрезок [a,b] разбит точками на n частичных отрезков imge016.gif (236 bytes). Сплайном степени m называется функция imge018.gif (230 bytes), обладающая следующими свойствами:

1) функция image019.gif (230 bytes) непрерывна на отрезке [a,b] вместе со своими производными до некоторого порядка  p.

2) на каждом частичном отрезке image021.gif (264 bytes) функция совпадает с некоторым алгебраическим многочленом image023.gif (240 bytes) степени m.

Разность m-p между степенью сплайна и наивысшим порядком непрерывной на отрезке [a,b] производной называют дефектом сплайна. Кусочно-линейная функция является сплайном первой степени с дефектом, равным единице. Действительно, на отрезке [a,b] сама функция image025.gif (225 bytes) (нулевая производная) непрерывна. В то же время на каждом частичном отрезке image026.gif (225 bytes) совпадает с некоторым многочленом первой степени.

Наиболее широкое распространение получили сплайны 3 степени (кубические сплайны) image028.gif (225 bytes) с дефектом равным 1 или 2. Система для осуществления сплайн-интерполяции кубическими полиномами предусматривает несколько встроенных функций. Одна из них рассмотрена в примере.

 

Погрешность приближения кубическими сплайнами. Пусть функция f имеет на отрезке [a,b] непрерывную производную четвертого порядка и image030.gif (597 bytes). Тогда для интерполяционного кубического сплайна справедлива оценка погрешности: image032.gif (739 bytes).

 

 

Наверх

10. Решение задачи Коши одношаговыми методами.

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

Пусть дано дифференциальное уравнение первого порядка image015.gif(302 bytes) . Требуется найти функцию image017.gif(228 bytes) , удовлетворяющую при image019.gif(220 bytes) дифференциальному уравнению и при image021.gif(222 bytes) начальному условию image023.gif(296 bytes)

Теорема существования и единственности задачи Коши. Пусть функция image038.gif(255 bytes) определена и непрерывна на множестве точек image040.gif(447 bytes) . Предположим также, что она удовлетворяет условию Липшица: image042.gif(566 bytes) для всех image044.gif(269 bytes) и произвольных image046.gif(193 bytes), image051.gif(196 bytes) , где image053.gif(179 bytes) - некоторая константа (постоянная Липшица). Тогда для каждого начального значения image055.gif(196 bytes) существует единственное решение image057.gif(228 bytes) задачи Коши, определенное на отрезке image059.gif(249 bytes) .

Геометрически задача интегрирования дифференциальных уравнений состоит в нахождении интегральных кривых, которые в каждой своей точке имеют заданное направление касательной. Заданием начального условия мы выделяем из семейства решений ту единственную кривую, которая проходит через фиксированную точку image061.gif(265 bytes)

 

Численное решение задачи Коши методом Эйлера.

Численное решение задачи Коши состоит в построении таблицы приближенных значений image063.gif(284 bytes)в точках image065.gif(265 bytes) . Точки image067.gif(277 bytes), image069.gif(268 bytes) называются узлами сетки, а величина image071.gif(182 bytes) - шагом сетки. В основе построения дискретной задачи Коши лежит тот или иной способ замены дифференциального уравнения его дискретным аналогом. Простейший метод основан на замене левой части уравнения правой разностной производной: image073.gif(453 bytes) . Разрешая уравнение относительно image075.gif(211 bytes) , получаем расчетную формулу метода Эйлера: image077.gif(394 bytes), image079.gif(268 bytes)

Численный метод называется  явным, если вычисление решения в следующей точке image082.gif(211 bytes)осуществляется по явной формуле. Метод называется одношаговым, если вычисление решения в следующей точке image084.gif(209 bytes) производится с использованием только одного предыдущего значения image086.gif(189 bytes) . Метод Эйлера является явным одношаговым методом. 

 

Оценка погрешности метода Эйлера.

Локальной погрешностью метода называется величина image089.gif(349 bytes) . Найдем величину локальной погрешности метода Эйлера: image091.gif(941 bytes) , при условии, что image093.gif(282 bytes) . Другими словами image095.gif(201 bytes) погрешность, которую допускает за один шаг метод, стартующий с точного решения. Глобальной погрешностью   (или просто погрешностью) численного метода называют сеточную функцию image097.gif(199 bytes) со значениями image099.gif(312 bytes) в узлах. В качестве меры абсолютной погрешности метода примем величину image101.gif(601 bytes) . Можно показать, что для явных одношаговых методов из того, что локальная погрешность имеет вид image103.gif(294 bytes) следует, что image105.gif(326 bytes) , где image107.gif(185 bytes) и M - некоторые константы. Таким образом, метод Эйлера является методом первого порядка точности. Для нахождения решения задачи Коши с заданной точностью image109.gif(177 bytes) требуется найти такое приближенное решение image111.gif(204 bytes) , для которого величина глобальной погрешности image113.gif(279 bytes) . Так как точное решение задачи неизвестно, погрешность оценивают с помощью правила Рунге.

Правило Рунге оценки погрешностей.  Для практической оценки погрешности проводят вычисления с шагами h и h/2. За оценку погрешности решения, полученного с шагом h/2, принимают величину, равную image115.gif(618 bytes) , где p - порядок метода.

 

Модификации метода Эйлера.

Метод Эйлера обладает медленной сходимостью, поэтому чаще применяют методы более высокого порядка точности. Второй порядок точности по image118.gif(182 bytes) имеет усовершенствованный метод Эйлера : image119.gif(702 bytes). Этот метод имеет простую геометрическую интерпретацию. Метод Эйлера называют методом  ломаных, так как интегральная кривая на отрезке image120.gif(262 bytes) заменяется ломаной с угловым коэффициентом image121.gif(278 bytes) . В усовершенствованном методе Эйлера интегральная кривая на отрезке image122.gif(263 bytes) заменяется ломаной с угловым коэффициентом, вычисленным в средней точке отрезка image123.gif(333 bytes) . Так как значение image124.gif(230 bytes) в этой точке неизвестно, для его нахождения используют метод Эйлера с шагом image125.gif(213 bytes)

 

Еще одна модификация метода Эйлера второго порядка - метод Эйлера-Коши:

image126.gif(829 bytes)

 

Решение систем дифференциальных уравнений методом Эйлера.

Пусть требуется решить нормальную систему дифференциальных уравнений: 

image127.gif(801 bytes)

с начальными условиями: image128.gif(310 bytes), image129.gif(318 bytes) , ..., image130.gif(319 bytes)

Эту систему в векторной форме можно записать в виде: image131.gif(313 bytes), image132.gif(309 bytes) . Здесь image133.gif(627 bytes), image134.gif(805 bytes), image135.gif(510 bytes) . Расчетная формула метода Эйлера имеет вид: image136.gif(410 bytes), image137.gif(268 bytes).

Теги

      21.12.2019

      Комментарии