Классификация и структура систем компьютерной математики

Министерство образования и науки РФ Государственное образовательное учреждение высшего профессионального образования ―Новгородский государственный университет имени Ярослава Мудрого» Институт электронных и информационных систем Кафедра проектирования и технологии радиоаппаратуры ЧИСЛЕННЫЕ МЕТОДЫ Конспект лекций Великий Новгород 2014 Лекция 1 КЛАССИФИКАЦИЯ И СТРУКТУРА СИСТЕМ КОМПЬЮТЕРНОЙ МАТЕМАТИКИ Содержание  Средства компьютерной математики.  Классификация компьютерных математических систем  Структура систем компьютерной математики  Системы компьютерной математики  Особенность численных методов по сравнению с аналитическими Введение С развитием средств вычислительной техники возникло новое научное и прикладное направление, получившее название компьютерной математики. Следует отметить, что термин «компьютерная математика» — это, прежде всего, собирательное название для целой серии математических систем, аккумулирующих многовековые знания человечества в области математических методов и расчетов, выполняемых как в численной, так и в аналитической и графической формах. Такие системы активно помогают массовому пользователю в его учебе или в проведении научной, педагогической и творческой деятельности. 1.1. Компоненты компьютерной математики Как и всякую информационную науку, компьютерную математику можно определить как совокупность теоретических, методических, аппаратных и программных средств. Эта область знаний обеспечивает эффективное выполнение с помощью компьютеров всех видов математических вычислений, как в автоматическом, так и диалоговом режимах с высокой степенью их визуализации. Аппаратные средства компьютерной математики используются в самых современных микрокалькуляторах, в математических сопроцессорах, графических процессорах и, наконец, в самих микропроцессорах, на которых строятся современные персональные компьютеры. Они включают в себя реализации новых команд МП, в том числе относящиеся к технологии MMX и средствам осуществления быстрой трехмерной графики. Для обычных пользователей такие средства практически закрыты, а их создание — прерогатива фирм, выпускающих МП. Теоретические средства включают в себя многие широко используемые математические методы:  интерполяцию и аппроксимацию функций;  численное решение трансцендентных и алгебраических уравнений;  численное интегрирование и дифференцирование;  обработку больших массивов данных; 2  дискретное преобразование Фурье;  двухмерную и трехмерную графику;  и т.д. Программные системы. Как всякое востребованное новое научное направление компьютерная математика бурно развивается. В ее недрах стали появляться присущие только ей отличительные продукты массового спроса — это, прежде всего, универсальные и специализированные программные системы компьютерной математики.  Система компьютерной алгебры MuPAD  Система MathCAD  Eureka  Derive  MATLAB (MathSoft Inc.)  Системы символьной математики  Mathematica  Maple V R4/R5 Большинство из перечисленных систем компьютерной математики имеет свое программное обеспечение в виде специализированного языка и пакетов расширений, позволяющих расширять функциональные возможности указанных систем. 1.2. Классификация компьютерных математических систем На сегодняшний день имеющиеся компьютерные математические системы можно (достаточно условно) разделить на 7 основных классов: 1. Системы для численных расчетов 2. Табличные процессоры 3. Матричные системы 4. Системы для статистических расчетов 5. Системы для специальных расчетов (САПР) 6. Системы для аналитических расчетов (компьютерной алгебры) 7. Универсальные системы Каждая из математических систем имеет определенные специфические свойства, которые предназначены для решения конкретных математических задач. Однако следует подчеркнуть, что отнесение современных компьютерных математических систем к одному классу практически невозможно. Хотя при создании системы ее предназначают для выполнения задач определенного класса, со временем в процессе модернизации ее характеристики приобретают универсальный характер. Поэтому можно говорить лишь о преимущественных сферах применения тех или иных математических систем. В 3 качестве примера, иллюстрирующего это тезис можно привести компьютерную систему Matlab. К сожалению, на нашем рынке и мировом рынках массовые системы компьютерной математики представлены только зарубежными программами. Это связано с тем, что современные программы этого класса относятся к числу наиболее сложных программных продуктов, требующих для своей разработки больших интеллектуальных, трудовых и финансовых затрат. Помимо указанного деления компьютерных систем на классы, правомерно также деление этих систем по сложности решаемых ими задач. Так системы Derive и MuPAD можно отнести к системам начального уровня, которые ориентированы на решение задач школьного образования и применения их студентами младших курсов вуза. К системам среднего класса можно отнести модернизированную систему MuPAD и ставшую весьма популярной систему MathCAD. Высший класс представлен системами Mathematica 2/3 и Maple V R4/R5. А такого «монстра» среди систем компьютерной математики, как матричную систему MATLAB 7 с ее многочисленными пакетами расширения и трудно укладываемой в нашем сознании стоимостью можно отнести к особо элитным и поэтому дорогим системам для избранных и весьма придирчивых пользователей. Это как бы «Мерседес 600» в мире математических систем. 1.3. Структура систем компьютерной математики Каждая система компьютерной математики может иметь нюансы в своей архитектуре и структуре. Тем не менее, большинство современных универсальных компьютерных систем базируется на следующей типовой структуре: Центральное место занимает ядро системы. Оно представляет множество заранее откомпилированных функций и процедур, представленных в машинных кодах и обеспечивающих набор встроенных функций и операторов системы. Это набор должен быть функционально полным. Например, в него должны входить как минимум все основные элементарные функции. Роль ядра особенно велика в системах символьной математики, где в ядре хранятся многие сотни, а то и тысячи правил преобразования математических выражений. Ядро математических систем тщательно оптимизируется поскольку от скорости его работы зависит скорость вычислений выполняемых данной системой 4 компьютерной математики. Этому способствует и предварительная компиляция ядра. Функции и процедуры включенные в ядро выполняются предельно быстро. С этой точки зрения в ядро было бы выгодно включать как можно больше вычислительных средств. Однако это приводит к замедлению поиска нужных средств из-за возрастания их числа, увеличению времени загрузки ядра и к другим нежелательным последствиям. Для оптимизации скорости работы объем ядра ограничивают, а расширение функциональных возможностей достигается за счет добавления к нему библиотек более редких процедур и функций, к которым обращается пользователь в случае их отсутствия в ядре. Некоторые системы допускают модернизацию библиотек и их расширение силами самих пользователей (m-файлы в MATLAB). Кардинальное расширение возможностей систем и их адаптация к реальным задачам, относящимся к разным областей знаний, достигается за счет пакетов расширения систем. Эти пакеты пишутся как правило на собственном языке программирования той или иной системы, что делает их подготовку обычными пользователями. Хотя как правило в базовую поставку включаются профессионально подготовленные фирменные пакеты расширения. Справочная система обеспечивает получение оперативных справок по вопросам работы с системами компьютерной математики, включая и примеры такой работы. Интерфейс дает пользователю возможность обращаться к ядру системы со своими запросами и получать результат решения на экране дисплея. Следует отметить, что большинство систем компьютерной математики имеют интерфейс, поддерживающий интерфейс ОС Windows. 1.4. Системы компьютерной математики 1.4.1. Встроенные калькуляторы Простейшими системами для численных расчетов являются встроенные в ОС микрокалькуляторы, которые не стоит путать с реальными калькуляторами, носимыми в кармане пиджака. Примером такого калькулятора является программное средство калькулятор, встроенное в ОС Windows, которое может работать в стандартном либо инженерном (научном) режимах. Может возникнуть вопрос, зачем работать с довольно примитивным калькулятором, когда под Windows разработан уже целый ряд достаточно мощных математических систем? Однако, как показывает опыт, именно при работе с такими системами нередко возникает нужда в подручном калькуляторе для проведения простейших вычислений, например по ходу ввода данных или обработке отдельных результатов. Примечание. Подобный режим простого калькулятора поддерживают и многие системы компьютерной математики, например Matlab. Работа с калькулятором отражает еще одну важную особенность Windows — обмен данными между отдельными приложениями. Это позволяет повысить производительность труда пользователя. Числа с экрана калькулятора можно скопировать в буфер и переслать, например в текстовый редактор WordPad или 5 Word. А можно поступить и наоборот. Таким образом, даже в таком калькуляторе реализована идея интеграции математических систем друг с другом. 1.4.2. Табличные процессоры В настоящее время известно множество программ для обработки табличных данных — табличных процессоров. Начало им было положено созданием электронных таблиц SuperCalc, которая применялась еще на старых 8-разрядных домашних компьютерах. К новейшим табличным процессорам принадлежат приложение Excel 2007, входящее в состав пакета офисных программ Office (последняя версия Office 2007) и приложение Calc, входящее в бесплатный пакет офисных программ OpenOffice.org 3.0 В настоящее время из-за своей специфики (необходимостью представления информации в табличной форме) эти пакеты в основном используются для финансово-экономических расчетов. 1.4.3. Системы для статистических расчетов Особую разновидность математических систем составляют программы, предназначенные для проведения статистических расчетов: StatGraphics Plus, Statistica, SPSS и другие. Есть среди этих программ и российская программа STADIA [х], созданная в МГУ. Интерфейс таких программ напоминает интерфейс программы Excel. Т.е. пользователю сначала необходимо ввести данные в виде электронной таблицы или загрузить их в таблицу с накопителей информации. Большинство вычислений выполняется по следующему алгоритму:  ввод данных в таблицу,  выделение нужных данных,  вызов команды, требующейся для выполнения необходимых вычислений (например, регрессии, корреляции, обработки временных рядов и т.д.). Главное отличие статистических программ от табличных процессоров заключается в большем числе встроенных специальных статистических функций, позволяющих выполнять без программирования огромное количество статистических вычислений, представляя их результаты в графической или иной форме. В повседневной жизни для выполнения статистических расчетов более удобно использование универсальных пакетов, которые менее перегружены специальными функциями, доступными только опытным специалистам по статистике. 1.4.4. Системы аналитических вычислений Системы аналитических вычислений (компьютерной или символьной алгебры) представляют собой новейшее направление развития систем компьютерной математики. Основное их достоинство заключается в возможности выполнения вычислений в аналитическом виде, а также в возможности проведения арифметических и многих иных вычислений практически с любой желаемой точностью и без ограничений по максимальным (минимальным) значениям чисел. 6 Системы символьной математики представляют собой наиболее актуальное и интересное направление компьютерной математики. Они уже сейчас делают то, что пару десятков лет тому назад казалось чистой фантастикой — выполняют сложнейший математические вычисления, в прошлом доступные только человеку. ЭВМ и программные системы, производящие символьные вычисления и способные выдавать результаты в виде аналитических формул, известны довольно давно. Лидирующая роль в разработке таких систем у нас принадлежала школе советского академика В.М. Глушкова, под руководством которого была создана серия инженерных ЭВМ серии Мир с языком Аналитик для проведения символьных вычислений. К сожалению, это направление не было поддержано в должной мере и лидерство перешло к зарубежным разработчикам таких средств. За рубежом был создан ряд таких систем, из которых наибольшее развитие в последнее время получили Maple V и Mathematica. Осознание важности компьютерной алгебры привело к тому, что ее средства в последствии были включены в более поздние версии наиболее серьезных систем для численных расчетов: Mathcad и MATLAB. Что превратило их в мощные и гибкие универсальные компьютерные системы. Отличительной чертой систем компьютерной алгебры является возможность вычисления математических выражений в общем виде. Например, если попытаться выполнить в обычной программе вычисление выражения типа sin2 (x) +cos2 (x) в общем виде для любого x — например, взять и поставить справа от выражения знак равенства, то при использовании обычных языков программирования или математических систем для численных расчетов будет выведено сообщение об ошибке — напоминание о том, что переменная x не определена. Это вполне естественно и совершенно тривиально. При использовании систем аналитических вычислений результатом будет, как и положено, единица. В этом принципиальное отличие данных систем от систем численного счета. Разумеется, это тривиальный пример. Куда важнее, что они способны вычислять аналитически производные и интегралы, выполнять подстановки одних сложных математических выражений в другие, выполнять математические преобразования (упрощение формул) и тому подобное. Словом выполнять «человеческую» работу. Численные методы решения математических задач при всей их огромной важности страдают двумя серьезными пороками: они дают только частные решения и всегда с конечной погрешностью. В противовес им символьные (аналитические) методы решения задач дают общие (формульные) результаты и без погрешности. Хотя, объективности ради. Следует сказать, что благодаря развитию и применению современных адаптивных и многовариантных численных методов вопрос о погрешности решений на практике уже решен и для численных методов. 1.4.5. Системы для специальных расчетов (САПР) Имеется большое число программ, изначально ориентированных на некоторые специализированные виды расчетов, связанных с особенностью их реализации для определенных предметных областей. 7 Из большого многообразия таких программ для нас интерес представляю два класса таких программ:  Программы для рисования графиков  САПР Из первого вида программ особо выделяется лучшая в этом классе — программа AXUM 5.0/6.0, созданная фирмой MathSoft Inc., и прекрасно интегрирующаяся с программой Mathcad. Эта программа обеспечивает около 70 типов графиков, многие из которых могут комбинироваться, богатейшее разнообразие вариантов графической визуализации результатов вычислений. К особому классу систем компьютерной математики относятся различные системы математического моделирования (САПР).  Например, для моделирования блочно заданных произвольных систем и устройств (системного моделирования) служит приложение Simulink, вошедшее в состав последних версий системы MATLAB.  Высокой степени совершенства достигли современные программы проектирования и моделирования электронных схем, среди которых можно выделить: MicroCAP, OrCAD, DesignLab, Electronics Workbench, AIM-Spice и другие. Они позволяют задавать в графическом виде электронные схемы, после чего автоматически формируют и решают весьма громоздкие системы алгебраических и дифференциальных уравнений. Результаты моделирования представляются в очень удобной графической форме в виде осциллограмм, спектрограмм, частотных характеристик и так далее. К сожалению, используемые в этих программах математические методы полностью скрыты от пользователя. Специфика использования этих программных продуктов, являющихся неотъемлемой частью арсенала компьютерной математики, будет рассмотрена позже в ряде специальных дисциплин, связанных с проектированием микросхем. 1.4.6. Универсальные компьютерные системы К этому классу можно отнести две системы:  Математическая система Mathcad.  Система MATLAB. Кратко перечислим особенности последней, поскольку данный курс базируется на этой системе. Система MATLAB Ранние версии систем MATLAB можно отнести к специализированным (матричным) системам, базирующимся на матричном аппарате. У этих систем единичное числовое значение воспринимается как элемент матрицы 11. Практически все функции системы (включая элементарные) определены как матричные — способные обрабатывать массивы. Название системы является сокращением слов MATrix LABoratory — матричная лаборатория. 8 В 90-х годах системы класса MATLAB интенсивно развивались в направлении универсализации, поэтому сейчас относить их только к матричным системам не совсем корректно. Последние версии этих систем (7.0 и 8.0) имеют не только самое большое количество матричных операторов, но и самую простую адаптацию к решению задач пользователя. Любое новое определение (расширение) системы задается в виде некоторого М-файла и в дальнейшем может использоваться наряду с встроенными в систему функциями без какого-либо объявления. Кроме того, для данной системы разработаны десятки пакетов расширения — от пакета символьной математики Symbolic до пакета имитационного моделирования блочно заданных систем Simulink. Так называемая дескрипторная графика этого пакета, основанная на принципах объектно-ориентированного программирования, добавляет мощные средства визуализации вычислений к огромному числу математических возможностей системы, среди которых много уникальных. 1.5. Особенность численных методов по сравнению с аналитическими Классические методы дифференцирования и интегрирования могут применяться для аналитически заданных функций. В тоже время на практике при проведении измерений или расчетах на компьютере результаты расчетов представляются в табличном либо графическом виде (так называемые сеточные функции), т.е. наборов дискретных значений аргументов. Дифференцирование и интегрирование сеточных функций требует разработки так называемых численных методов, учитывающих специфику таких функций. Прогресс в области создания изделий микроэлектроники оказал значительное влияние в первую очередь на методы их проектирования. Создание первых ИС (в монолитном и гибридном исполнении) поставило под сомнение традиционные методы их проектирования. Известно, что ручные методы проектирования (под ними будем подразумевать анализ упрощенных аналитических выражений, описывающих поведение электрической или электронной схемы в различных режимах) использовались в сочетании с макетированием или физическим моделированием, причем электрические или электронные схемы рассматривались как системы с сосредоточенными параметрами. Эволюция технологии ИС, сопровождающаяся непрерывным снижением величины топологической нормы, привела к тому, что даже достаточно точные модели первых типов интегральным схем оказались неработоспособными для новых поколений ИС. Благодаря этому, а отчасти из-за требования к постоянному сокращению сроков проектирования (из-за морального старения) при росте сложности ИС, приводит к тому, что ручные методы оказываются экономически крайне неэффективными. Поэтому возникло новое направление, занимающееся развитием и внедрением методов проектирования схем с помощью ЭВМ. Основу их составляет разработка точных математических моделей, представляющих собой систему математических выражений вида F= (x, dx/dt, B) 9 или в более привычной форме: I = (V, dV/dt, B). Это направление получило название автоматизированного проектирования. Сравним алгоритмы аналитических методов с численными на примере решения системы линейных уравнений вида Ax=b, известными вам еще со школы. Наиболее простой способ — использование обратной матрицы A-1 A-1Ax= A-1b, что приводит к решению в виде x= A-1b, Другой распространенный способ — метод исключения Гаусса. Он позволяет сократить объем вычислений по сравнению с предыдущим методом в два раза и выполняется в два этапа. Сначала матрица коэффициентов модифицируется в форму верхней треугольной матрицы (матрицы с нулевыми элементами ниже главной диагонали). Второй этап, известный как метод обратной подстановки, включает в себя решение последнего (n-го) уравнения для определения xn. Затем последовательно определяются xn-1, xn-2 и т.д. вплоть до x1. Для решения системы nго порядка методом Гаусса требуется выполнить следующее количество операций: n3/3 + n2 - n/3. Именно таким способом оценивается эффективность используемого численного метода для решения конкретной задачи. Сложность анализируемых схем постоянно возрастает, что приводит к увеличению числа решаемых уравнений и значительному росту машинного времени на проведение однократного анализа схем (расчета рабочей точки). КОНТРОЛЬНЫЕ ВОПРОСЫ 1. 2. 3. 4. 5. Как вы понимаете термин компьютерная математика? На каких компонентах базируется компьютерная математика? Нарисуйте базовую структурную схему систем компьютерной математики. Дайте классификацию систем компьютерной математики. между численными и аналитическими методами решения математических задач. 10 Лекция 2 ПРИБЛИЖЕННЫЕ ЧИСЛА, РАБОТА С НИМИ И ПОГРЕШНОСТИ Содержание         Приближенные числа Правила записи приближенных чисел Округление приближенных чисел. Абсолютная и относительная погрешность вычисления функции одной переменной. Абсолютная и относительная погрешность вычисления функции нескольких переменных Абсолютная и относительная погрешность вычисления суммы и разности приближенных чисел Абсолютная и относительная погрешность вычисления произведения и частного приближенных чисел Источники погрешности решения задачи на ЭВМ 2.1. Приближенные числа Пусть X – некоторая величина, истинное значение которой известно или неизвестно и равно x*. Число x, которое можно принять за значение величины X, мы будем называть ее приближенным значением или просто приближенным числом. Число x называют приближенным значением по недостатку, если оно меньше истинного значения (x < x* ), и по избытку, если оно больше (x > x* ). Например, число 3,14 является приближенным значением числа π по недостатку, а 2,72 – приближенным значением числа е по избытку. Абсолютная погрешность приближенного числа есть абсолютная величина разности между истинным значением величины и данным ее приближенным значением. Δx = |x* – x| Поскольку истинное значение величины обычно остается неизвестным, неизвестной остается также и абсолютная погрешность. Вместо нее приходится рассматривать оценку абсолютной погрешности, так называемою предельную абсолютную погрешность, которая означает число, не меньшее абсолютной погрешности (далее, в том случае, если это не принципиально, будем под абсолютной погрешностью понимать именно предельную абсолютную погрешность). Абсолютная погрешность приближенного числа не в полной мере характеризует его точность. Действительно, погрешность в 0,1 г слишком велика при взвешивании реактивов для проведения микро-синтеза, допустима при взвешивании 100 г колбасы, и не может быть замечена при измерении массы, например, железнодорожного вагона. Более информативным показателем точности приближенного числа является его относительная погрешность. 11 Относительной погрешностью δx приближенного значения величины X называют абсолютную величину отношения его абсолютной погрешности к истинному значению этой величины. Часто относительную погрешность выражают в процентах. положительности абсолютной погрешности можно записать: C учетом δx =Δx / |x* | Ввиду того, что фактически вместо абсолютной погрешности приходится рассматривать предельную, относительную погрешность также заменяют предельной относительной погрешностью, которая означает число, не меньшее относительной погрешности. Более того, при отыскании предельной относительной погрешности приходится заменять неизвестное истинное значение величины x* приближенным – x. Последняя замена обычно не отражается на величине относительной погрешности ввиду близости этих значений и малости абсолютной погрешности. δx =Δx / |x | Например, для приближенного значения π = 3,14 предельная абсолютная погрешность составляет 0,0016, а относительная – 0,00051 или 0,051%. Выражение относительной погрешности в процентах иногда называют процентной погрешностью. 2.2. Правила записи приближенных чисел Приближенные числа принято записывать таким образом, чтобы вид числа показывал его абсолютную погрешность, которая не должна превосходить половины единицы последнего разряда, сохраняемого при записи. Например, запись 3,1416 означает, что абсолютная погрешность этого приближенного числа не превосходит 0,00005. Для числа 370 абсолютная погрешность не превосходит 0,5. Если это число имеет большую точность, например если абсолютная погрешность меньше 0,05, то следует писать уже не 370, а 370,0. Таким образом, приближенные числа 37•101; 370; 370,0; 370,00 имеют различную степень точности: их предельные абсолютные погрешности составляют соответственно 5; 0,5; 0,05 и 0,005. Для больших чисел абсолютные погрешности могут иметь порядок единиц, десятков, сотен и т. п. Сохраняя и для таких чисел упомянутое выше правило, не следует выписывать все его цифры. Так, число двести семьдесят пять тысяч с абсолютной погрешностью, не превосходящей 500, надо писать в виде 275•103. Запись вида 275 000 применять нельзя, ибо она означала бы предельную абсолютную погрешность 0,5. Если абсолютная погрешность величины Х не превышает одной единицы разряда последней цифры приближенного числа х, то говорят, что у числа х все знаки верные. Приближенные числа следует записывать, сохраняя только верные знаки. Если, например, абсолютная погрешность числа 12400 равна 100, то это число должно быть записано, например, в виде 124•102 или 0,124•105. Оценить погрешность приближенного числа можно, указав, сколько верных значащих цифр оно содержит. При подсчете значащих цифр не считаются нули с левой стороны числа. 12 Например, в числе 0,0645 – 3 значащие цифры, а у числа 34,5400 – 6 значащих цифр. Высказанные правила записи приближенных чисел применяются для записи исходных данных, полученных в результате измерений, и в математических таблицах. Абсолютная погрешность при этом не выписывается, и всегда предполагается, что она не превосходит половины единицы последнего разряда, сохраняемого при записи. В указанной форме записи все цифры, таким образом, являются верными. В окончательных результатах расчета принято записывать числа с одной сомнительной цифрой (т. е. сохраняя следующую за верной). При этом следует указывать предельную абсолютную погрешность, выписывая ее с одной значащей цифрой. Для этого погрешность округления числа прибавляют к предельной абсолютной погрешности и результат округляют в сторону увеличения ). 2.3. Округление приближенных чисел Если приближенное число содержит лишние (или неверные) знаки, то его следует округлить, отбрасывая излишние цифры и руководствуясь следующим известным правилом округления: если первая из отброшенных цифр 4 или меньше, то последняя оставшаяся цифра сохраняется без изменения; если первая из отброшенных цифр 5 или больше, то последняя оставшаяся цифра увеличивается на единицу. Исключением из этого правила является случай, когда отбрасывается только пятерка или же пятерка с нулями. Здесь принято сохранять последнюю оставшуюся цифру без изменения, если она четная, и увеличивать ее на единицу до четной, если она была нечетная. При проведении вычислений рекомендуется сохранять максимальное число значащих цифр для всех промежуточных результатов. Округляется только конечный результат в соответствии с оцененной предельной абсолютной погрешностью. Если погрешность результата не оценивается, то можно руководствоваться следующими общими правилами округления результатов действий:  При сложении и вычитании приближѐнных чисел в результате следует сохранять столько десятичных знаков, сколько их в приближѐнном данном с наименьшим числом десятичных знаков.  При умножении и делении в результате следует сохранять столько значащих цифр, сколько их имеет приближѐнное данное с наименьшим числом значащих цифр.  При возведении в квадрат или куб в результате следует сохранять столько значащих цифр, сколько их имеет возводимое в степень приближѐнное число (последняя цифра квадрата и особенно куба при этом менее надежна, чем последняя цифра основания ).  При извлечении квадратного и кубического корней в результате следует брать столько значащих цифр, сколько их имеет приближѐнное значение 13 подкоренного числа (последняя цифра квадратного и особенно кубического корня при этом более надѐжна, чем последняя цифра подкоренного числа). 2.4. Абсолютная и относительная погрешность вычисления функции одной переменной Важной проблемой при проведении вычислений с использованием приближенных чисел является вопрос о влиянии погрешности исходных данных на погрешность полученного результата. Установим основные свойства распространения абсолютной и относительной погрешности. Для начала рассмотрим случай вычисления функции одного аргумента, который является приближенным числом. Решение этих задач дается следующей основной теоремой. Теорема. Предельная абсолютная погрешность вычисления функции равна произведению абсолютной величины ее производной на предельную абсолютную погрешность аргумента. Доказательство. Пусть число x является приближенным значением величины X с абсолютной погрешностью Δx: , или можно записать , если считать, что величина Δx является знаковой, то есть, может быть как положительным, так и отрицательным числом, потому что x может быть приближением как по недостатку, так и по избытку. Обозначим абсолютную погрешность функции через Δy . Ввиду малости Δx мы можем заменить приращение функции ее дифференциалом. Тогда получим , откуда , что и требовалось доказать. Рассмотрим относительную погрешность вычисления функции одной переменной. Обозначим предельную относительную погрешность аргумента через δx, а функции – δy, тогда, учитывая, что , 14 получим . Так как по правилам дифференцирования сложной функции , то полученное выражение можно записать так: . В таком виде выражение удобно для приложения к легко логарифмируемым функциям. Рассмотрим, например, погрешность вычисления степенной функции: . Здесь , т. е. предельная относительная погрешность степени равна предельной относительной погрешности основания, умноженной на абсолютную величину показателя степени. 2.5. Абсолютная и относительная погрешность вычисления функции нескольких переменных Полученные результаты легко обобщаются на случай функций нескольких аргументов . Абсолютная погрешность результата вычисления функции нескольких приближенных чисел равна сумме произведений модуля частной производной функции на абсолютную погрешность приближенного числа. 15 Для относительной погрешности вычисления функции нескольких приближенных чисел получим выражение, так же аналогичное случаю функции одной переменной: . 2.6. Абсолютная и относительная погрешность вычисления суммы и разности приближенных чисел Рассмотрим несколько примеров функций частного вида: , таким образом, , т. е. предельная абсолютная погрешность как суммы, так и разности нескольких приближенных чисел равна сумме предельных абсолютных погрешностей слагаемых. При определении относительной погрешности суммы и разности результаты будут различны. Для суммы n положительных чисел относительная погрешность будет равна . Как видно, величина относительной погрешности в этом случае зависит не только от погрешностей слагаемых, но и от величины приближенных чисел. Можно воспользоваться следующими приемом для оценки границ относительной погрешности суммы. Выберем из всех относительных погрешностей наибольшую , тогда можно записать 16 Таким образом, предельная относительная погрешность суммы положительных чисел не превосходит максимальной относительной погрешности слагаемых, т.е. является оценкой сверху относительной погрешности суммы приближенных чисел. Аналогично, если выбрать наименьшую из относительных погрешностей слагаемых, можно показать, что она будет являться оценкой снизу относительной погрешности суммы приближенных чисел. Рассмотрим далее относительную погрешность вычисления разности двух приближенных чисел . Как видно, из полученного результата, относительная погрешность разности может очень существенно (в пределе неограниченно) возрастать при вычитании двух близких по значению чисел . Пример 2.1 Рассмотрим вычисление массы образца как разности масс сосуда с образцом (x1) и пустого сосуда (x2) x1 = 1,245 x1 = 0,0005 x2 = 1,235 x1 = 0,0005 Относительные погрешности примерно одинаковы и равны 0,04%. Найдем массу навески как разность двух чисел y = 1,245 - 1,235 = 0,010 . Абсолютная погрешность y = 0,0005 + 0,0005 = 0,001 . Таким образом, относительная погрешность результата составит y = | 0,001/ 0,01 | = 0,1 = 10% , что в 250 (!!!!) раз выше относительной погрешности исходных чисел. 2.7. Абсолютная и относительная погрешность вычисления произведения и частного приближенных чисел Рассмотрим произведение нескольких приближенных чисел. y = x1 · x2 ·…· xn . 17 Воспользуемся общей формулой для относительной погрешности, в логарифмическом виде, полученной ранее. Для этого найдем натуральный логарифм функции ln y = ln x1 + ln x2 +…+ xn и, затем, частные производные логарифма функции по всем xi. Для всех i имеем . Поэтому, окончательно, получим: . То есть, предельная относительная погрешность произведения приближенных чисел равна сумме предельных относительных погрешностей сомножителей. Для частного двух чисел вычисления будут аналогичными. y = x1 / x2 . Находим натуральный логарифм и его частные производные ln y = ln x1 - ln x2 , , . Следовательно, предельная относительная погрешность частного равна сумме предельных относительных погрешностей делимого и делителя. Пример 2.2. При определении силы сопротивления квадратной пластинки, поставленной перпендикулярно к воздушному потоку, используют формулу , где Ρ — сила сопротивления, S — площадь пластинки, ν — скорость воздушного потока и k — коэффициент пропорциональности. Известно, что величина k определена с относительной погрешностью 5%, S и v — с относительными погрешностями 1%. Требуется определить относительную погрешность величины Р. Решение Согласно выражению для предельной относительной погрешности произведения имеем , 18 откуда . Итак, для оценки погрешности мы получили следующие простые правила:  При сложении и вычитании абсолютные погрешности складываются;  При умножении и делении относительные погрешности складываются;  При возведении в степень относительные погрешности умножаются на абсолютную величину показателя степени;  При отыскании значения функции абсолютная погрешность функции равна произведению абсолютной погрешности аргумента на абсолютную величину производной Заметим еще, что при вычислении значений функции абсолютная погрешность может существенным образом зависеть от того, каким образом записана расчетная формула и какова последовательность операций в этой формуле. Формулы для вычислений надо стараться преобразовывать к такому виду, чтобы в них не было вычитания близких величин; последнее может привести к большой потере точности и большим относительным ошибкам. 2.8. Источники погрешности решения задачи на ЭВМ Численные методы решения различных нелинейных математических задач по своей природе являются приближенными в отличие от прямых методов, дающих точное решение. Поэтому с точки зрения точности результата использование прямых методов может показаться более предпочтительным. Однако на самом деле при решении задачи на компьютере ответ все равно будет содержать погрешность. В качестве основных источников погрешности обычно рассматривают три вида ошибок:  ошибки усечения,  ошибки округления,  ошибки распространения. Рассмотрим их подробнее. 2.8.1. Ошибки усечения Этот вид ошибок связан с погрешностью, заложенной в самой задаче. Он может быть обусловлен неточностью определения исходных данных. Например, если в условии задачи заданы какие-либо размеры, то на практике для реальных объектов эти размеры известны всегда с некоторой точностью. То же самое касается любых других физических параметров. Сюда же можно отнести неточность расчетных формул и входящих в них числовых коэффициентов. Большое число расчетных формул являются эмпирическими и дают результат с некоторой погрешностью, содержат подгоночные коэффициенты, обеспечивающие приемлемую ошибку в ограниченном диапазоне входных параметров. Поэтому, как 19 правило, если исходные данные известны с некоторой погрешностью, вряд ли стоит пытаться получить результат с меньшей погрешностью. 2.8.2. Ошибки распространения Данный вид ошибок связан с применением того или иного способа решения задачи. В ходе вычислений неизбежно происходит накопление или, иначе говоря, распространение ошибки. Помимо того, что сами исходные данные не являются точными, новая погрешность возникает при их перемножении, сложении и т. п. Накопление ошибки зависит от характера и количества арифметических действий, используемых в расчете. Обычно для решения одной и той же задачи может быть использован ряд различных методов решения. Например, систему линейных алгебраических уравнений можно решить методом Гаусса или через определители (методом Крамера). Теоретически оба метода позволяют получить точное решение. Однако на практике при решении больших систем уравнений метод Гаусса обеспечивает меньшую погрешность, чем метод Крамера, так как использует меньший объем вычислений. 2.8.3. Ошибки округления Это тип ошибок связан с тем, что истинное значение числа не всегда точно сохраняется компьютером. При сохранении вещественного числа в памяти компьютера оно записывается в виде мантиссы и порядка, примерно так же, как отображается число на дисплее калькулятора (см. рис. 2.1). Рис. 2.1. Структура записи вещественного числа Здесь R1, R2, R3, … Rn – разряды мантиссы, D1, D2, …, Dm – разряды порядка. На самом деле конечно, в отличие от дисплея калькулятора, мантисса и порядок числа, включая их знаки, в памяти компьютера хранятся в двоичном виде. Но для обсуждения природы ошибок округления это различие не столь принципиально. Понятно, что иррациональные числа такие, как π = 3,14159… и e = 2,712… не могут быть представлены в памяти компьютера в принципе. Однако же и рациональные числа, если количество их значащих цифр превышает число отведенных разрядов мантиссы (см. рис. 2.1), будут представлены не точно. При этом цифра последнего сохраняемого в ЭВМ разряда может быть записана с округлением или без него. Фактически при заданной структуре хранения числа компьютер может использовать не бесконечное, а конечное число рациональных чисел, которые вписываются в приведенную на рис. 2.1 схему. Поэтому любой входной параметр решаемой задачи, ее промежуточный результат и окончательной ответ всегда округляются до разрешенных в компьютере чисел. Следующий важный вывод касается диапазона представления чисел в ЭВМ. Если проводить рассуждения для десятичной системы счисления, то максимальное по модулю число, которое может быть представлено в соответствии со схемой на рис. 2.1, равно 20 ±X∞ = ± 999 … 9 × 10 +99 … 9 . Все числа, превышающие по модулю X∞, не представимы в ЭВМ и рассматриваются как машинная бесконечность. Если в ходе расчетов будет получен результат, превышающий X∞, то произойдет аварийное завершение вычислений по переполнению. Нетрудно убедиться опытным путем, что, например, в MATLAB верхний диапазон чисел соответствует X∞ ~ 10307. Минимальное по модулю число, сохраняемое в памяти компьютера, по схеме на рис. 2.1 равно ±X0 = ± 000 … 1 × 10 –99 … 9 . Числа, модуль которых меньше X0, воспринимаются ЭВМ как нуль, точнее как машинный нуль. Если при выполнении расчетов будет получен результат меньше, чем X0, то это будет воспринято как потеря порядка. Обычно в подобной ситуации результат полагается равным нулю, и вычисления продолжаются. На рис. 2.2 показана "машинная" числовая ось, на которой отмечены X0 и X∞ . Числа располагаются на оси неравномерно. Их плотность возрастает по мере приближения к нулю. Рис. 2.2. «Машинная» числовая ось На рис. 2.2 вблизи единицы отмечена небольшая область εм, которую называют машинное эпсилон. Параметр εм весьма важен, так как он характеризует относительную точность представления чисел в компьютере. В зависимости от способа округления чисел в ЭВМ величина εм определяется первым отбрасываемым или последним сохраняемым разрядом мантиссы. Следует иметь ввиду, что длина мантиссы в памяти компьютера устанавливается программно. Например, при выполнении расчетов на языке ФОРТРАН с использованием "обычной" точности (двоичная запись числа длиной четыре байта) εм ~ 10–7. При удвоенной длине мантиссы εм ~ 10–16. КОНТРОЛЬНЫЕ ВОПРОСЫ Вопрос 1. Объясните правила записи приближенных чисел. Вопрос 2. Сформулируйте правила округление приближенных чисел. Вопрос 3. Объясните смысл абсолютной и относительной погрешности и способы их вычисления для функции одной переменной. Вопрос 4. Объясните смысл абсолютной и относительной погрешности и способы их вычисления для функции нескольких переменных. Вопрос 5. Сформулируйте правила вычисления абсолютной и относительной погрешностей для суммы и разности приближенных чисел Вопрос 6. Сформулируйте правила вычисления абсолютной и относительной погрешностей для произведения и частного приближенных чисел. Вопрос 7. Сформулируйте источники погрешности решения задачи на ЭВМ. 21 Лекция 3 ЧИСЛЕННЫЕ МЕТОДЫ РЕШЕНИЯ НЕЛИНЕЙНЫХ УРАВНЕНИЙ Содержание  Определение простого и кратного корней.  Постановка задачи приближенного решения нелинейного уравнения.  Алгоритм метода дихотомии (деления отрезка пополам)  Теорема о сходимости метода хорд.  Алгоритм метода Ньютона.  Теорема о сходимости метода Ньютона  Алгоритм метода простой итерации  Теорема о сходимости метода простой итерации  Приведение уравнения к виду, удобному для итераций. 3.1. Введение Пусть рассматривается уравнение значение , при котором . Корнем уравнения называется . Корень называется простым, если , в противном случае корень называется кратным. Целое число m называется кратностью корня . Постановка задачи. точностью , если Вычисление для k=1,2,3-,m-1 и приближенного предполагает поиск найти такое значения значения корня с , что . Решение задачи разбивается на два этапа:  на первом этапе осуществляют локализацию корней,  на втором этапе производят итерационное уточнение корней. На этапе локализации корней находят достаточно узкие отрезки (или отрезок, если корень единственный), которые содержат один и только один корень уравнения f(x)=0. На втором этапе вычисляют приближенное 22 значение корня с заданной точностью. Часто вместо отрезка локализации достаточно указать начальное приближение к корню. Примечание. Для практического знакомства с процедурами локализации корней выполните лаб. работу № 2. 3.2. Метод дихотомии (деления отрезка пополам) Пусть [a,b] – отрезок локализации. Предположим, что функция f(x) непрерывна на [a,b] и на концах принимает значения разных знаков, т.е. f(a)* f(b) 0. Алгоритм метода дихотомии состоит в построении последовательности вложенных отрезков, на концах которых функция принимает значения разных знаков. Каждый последующий отрезок получают делением пополам предыдущего. Опишем один шаг итераций метода. Пусть на k-ом шаге найден отрезок такой, что отрезка . Если . Найдем середину , то - корень и задача решена. Если нет, то из двух половин отрезка выбираем ту, на концах которой функция имеет противоположные знаки: , если , , , если Критерий окончания итерационного процесса: если длина отрезка локализации меньше 2ε, то итерации прекращают и в качестве значения корня с заданной точностью принимают середину отрезка. 3.2.1. Теорема о сходимости метода дихотомии Пусть функция f(x) непрерывна на [a,b] и на концах принимает значения разных знаков f(a)* f(b) 0. Тогда метод сходится и справедлива оценка погрешности: 23 Примечание. Для практического знакомства с технологий численного решения нелинейного уравнения методом дихотомии выполните лаб. работу № 4. 3.3. Метод Ньютона (метод касательных) Расчетная формула метода Ньютона имеет вид: . Геометрически метод Ньютона означает, что следующее приближение к корню x(n+1) есть точка пересечения с осью ОХ касательной, проведенной к графику функции y = f(x) в точке (x(n), f(x(n))). 3.3.1. Теорема о сходимости метода Ньютона Пусть - простой корень уравнения f(x)=0, в некоторой окрестности которого функция дважды непрерывно дифференцируема. Тогда найдется такая малая ζ - окрестность корня начального приближения x(0) из , что при произвольном выборе этой окрестности итерационная последовательность метода Ньютона не выходит за пределы окрестности и справедлива оценка , где Критерий окончания итерационного , . процесса : при заданной точности  > 0 вычисления следует вести до тех пор пока не окажется выполненным неравенство . Примечание. Для практического знакомства с технологий численного решения нелинейного уравнения методом Ньютона выполните лаб. работу № 3. 24 Как указано в теореме, метод Ньютона обладает локальной сходимостью, то есть областью его сходимости является малая окрестность корня . Неудачный выбор может дать расходящуюся итерационную последовательность. 3.4. Метод простой итерации (метод последовательных повторений) Для применения метода простой итерации следует исходное уравнение f(x) = 0 преобразовать к виду, удобному для итерации x = (x). Это преобразование можно выполнить различными способами. Функция (x) называется итерационной функцией. Расчетная формула метода простой итерации имеет вид: x(n+1) = (x(n)). Теорема о сходимости метода простой итерации. Пусть в некоторой  - окрестности корня функция (x) дифференцируема и удовлетворяет неравенству |’(x)| q, где 0 q 1 — постоянная. Тогда независимо от выбора начального приближения из указанной  - окрестности итерационная последовательность не выходит из этой окрестности, метод сходится со скоростью геометрической последовательности и справедлива оценка погрешности: , . Критерий окончания итерационного процесса . При заданной точности ε >0 вычисления следует вести до тех пор, пока не окажется выполненным неравенство . 25 Если величина 0  q  0,5, то можно использовать более простой критерий окончания итераций: . Внимание. Ключевым моментом в применении метода простой итерации является эквивалентное преобразование уравнения к новой форме. Способ, при котором выполнено условие сходимости метода простой итерации, состоит в следующем: исходное уравнение приводится к виду x = x -  f(x). Предположим дополнительно, что производная f’ непрерывна и положительна на заданном отрезке, т.е. справедливо выполнение неравенства m f’(x) M на отрезке [a,b]. Тогда при выборе итерационного параметра  = 2/(m+M) метод сходится и значение еq  M m  1. M m Если же известна для производной только оценка сверху, то положим  = 1/ M и q = 1 – m/M – тоже менее единицы. Примечание. Для практического знакомства с технологий численного решения нелинейного уравнения методом простой итерации выполните лаб. работу № 3. Задачи для самостоятельного решения Задача 1. Определить количество корней уравнения и для каждого корня найти отрезки локализации: x4 – 1 – cos(x) = 0. Задача 2. Методом деления отрезка пополам найти корень уравнения x3 + 2x – 6 = 0 с точностью  = 0.3. Сколько нужно сделать итераций для получения точности  = 0.01? Задача 3. Методом простой итерации найти корни уравнения x4 – 4x – 4 = 0 c точностью  = 0.1 Указание: отрицательный корень найти простым преобразованием уравнения, а положительный корень найти методом простой итерации с оптимальным выбором итерационного параметра. 26 Задача 4. Записать расчетную формулу метода Ньютона и указать критерий окончания итераций для решения уравнения 10x + x2 -2 = 0. Вычислить два первых приближения к корню. Задача 5. Указать критерий окончания в методе Ньютона для решения задачи: методом Ньютона найти корень уравнения с k верными значащими цифрами. Контрольные вопросы Вопрос 1. Сформулируйте постановку задачи приближенного решения нелинейного уравнения и основные этапы ее решения. Вопрос 2. Докажите оценку погрешности метода дихотомии. Вопрос 3. Запишите расчетную формулу метода Ньютона и дайте геометрическую интерпретацию метода. Вопрос 4. Что такое итерационная функция? Вопрос 5. Выведете критерий окончания итераций для метода простой итерации из оценки погрешности. Вопрос 6. Можно ли найти кратный корень с помощью метода дихотомии? 27 Лекция 4 РЕШЕНИЕ СЛАУ ПРЯМЫМИ МЕТОДАМИ Содержание  Прямые и итерационные методы.  Метод Гаусса (схема частичного выбора).  Метод Холецкого.  Метод прогонки. 4.1. Прямые и итерационные методы Метод решения задачи называют прямым, если он позволяет получить решение после выполнения конечного числа элементарных операций. Метод решения задачи называют итерационным, если в результате получают бесконечную последовательность приближений к решению. Если эта последовательность сходится к решению задачи, то говорят, что итерационный процесс сходится. Эти методы будут рассмотрены в следующей лекции. К прямым методам решения относятся метод Гаусса и его модификации, метод Холецкого и метод прогонки. 4.2. Метод Гаусса с выбором главного элемента по столбцу (схема частичного выбора) В методе Гаусса для вычисления масштабирующих множителей требуется делить на ведущие элементы каждого шага. Если элемент равен нулю или близок к нулю, то возможен неконтролируемый рост погрешности. Поэтому часто применяют модификации метода Гаусса, обладающие лучшими вычислительными свойствами. В качестве примера рассмотрим одну из таких модификаций — метод Гаусса с выбором главного элемента по столбцу (схема частичного выбора) Здесь на k-ом шаге прямого хода в качестве ведущего элемента выбирают максимальный по модулю коэффициент при неизвестной xk в уравнениях с номерами i = k+1, ... , m. Затем уравнение, соответствующее выбранному коэффициенту с номером ik, меняют местами с к-ым уравнением системы для того, чтобы главный элемент занял место коэффициента akk(k-1). После этой перестановки исключение проводят как в схеме единственного деления. В этом случае все масштабирующие множители по модулю меньше единицы и схема обладает вычислительной устойчивостью. 28 Пример 1. Решение системы методом Гаусса с выбором главного элемента по столбцу. Пусть Ax = b, где A= , b= Прямой ход. Шаг 1. Максимальный по модулю элемент 1-го столбца a33=-27. Переставим 1-ое и 3 - е уравнения местами: A= , b= Вычислим масштабирующие множители 1 шага: и выполним преобразование матрицы и вектора: A1= b1= Шаг 2. Вычислим масштабирующие множители 2 шага: . Второй шаг не изменяет матриц: A2=A1, b2= b1. Шаг 3. Максимальный по модулю элемент 3 столбца a43= 43/3. Переставим 3 и 4 уравнения местами. A2 = , b2 = 29 Вычислим масштабирующие множители 3 шага: и выполним преобразование матрицы и вектора: A3 = , b3 = Обратный ход. Из последнего уравнения находим: x4 = 0. Из третьего уравнения системы находим . Из второго уравнения находим . Неизвестное x1 находим из первого уравнения: Ответ: . 4.3. Метод Холецкого Если матрица системы является симметричной и положительно определенной, то для решения системы применяют метод Холецкого (метод квадратных корней). В основе метода лежит алгоритм специального LU-разложения матрицы A, в результате чего она приводится к виду A = LLT . Если разложение получено, то как и в методе LU-разложения, решение системы сводится к последовательному решению двух систем с треугольными матрицами: Ly = b и LTx = y. Для нахождения коэффициентов матрицы L неизвестные коэффициенты матрицы LLT приравнивают соответствующим элементам матрицы A. Затем последовательно находят требуемые коэффициенты по формулам: , i = 2, 3, ..., m, , i = 3, 4, ..., m, 30 i = k+1, ... , m. Пример 2. Решение системы методом Холецкого Пусть A= , b= Находим элементы матрицы L: Таким образом разложение матрицы A имеет вид: Последовательно решаем системы Ly = b и LTx = y. Решением 1-ой системы является вектор , а решение 2-ой системы вектор .. Ответ: x1 = 6, x2 = -5, x3 = -4 31 4.4. Метод прогонки Если матрица системы является разреженной, то есть содержит большое число нулевых элементов, то применяют еще одну модификацию метода Гаусса - метод прогонки. Рассмотрим систему уравнений с трехдиагональной матрицей: Преобразуем первое уравнение системы к виду x1 = a1x2 +  1, где a1 = c1 /b1, 1 = d1 /b1. Подставим полученное выражение во второе уравнение системы и преобразуем его к виду x2 = a2x3 +  2 и т.д. На i-ом шаге уравнение преобразуется к виду xi = ai xi+1 +  i , где ai = -c1 /(bi + aiai-1),  i = (di - ai i-1) /(bi + aiai-1). На m-ом шаге подстановка в последнее уравнение выражения xm-1 = am-1 xm +  m-1 , дает возможность определить значение xm: xm =  m = (dm – am m-1) /(bm + amam-1). Значения остальных неизвестных находятся по формулам: xi = ai xi+1 +  i , i = m-1, m-2, ..., 1. Пример 4. Решение системы уравнений методом прогонки Прямой ход прогонки. Вычислим прогоночные коэффициенты: , , 32 , , , , Обратный ход прогонки. Находим значения неизвестных: , , Ответ: , . ЗАДАЧИ ДЛЯ САМОСТОЯТЕЛЬНОГО РЕШЕНИЯ Задача 1. Методом Гаусса с выбором главного элемента найти решение системы уравнений: Задача 2. Найти разложение матрицы: Задача 3. Можно ли применять метод Холецкого к системе уравнений, матрица которой имеет вид: Задача 4. Подсчитать количество арифметических действий в методе прогонки. КОНТРОЛЬНЫЕ ВОПРОСЫ Вопрос 1. Что такое прямые и итерационные методы. Вопрос 2 С какой целью применяют модификацию метода Гаусса - схему частичного выбора. Вопрос 3. Для каких систем уравнений применяют метод Холецкого. 33 Вопрос 4 Запишите формулы для нахождения решения после приведения системы к виду. Вопрос 5. Сформулируйте алгоритм метода прогонки. 34 Лекция 5 РЕШЕНИЕ СЛАУ ИТЕРАЦИОННЫМИ МЕТОДАМИ Содержание  Приведение системы к виду удобному для итераций  Достаточные условия сходимости  Критерий окончания итераций.  Метод Якоби  Метод Зейделя  Метод простой итерации. 5.1. Алгоритм и условие сходимости итерационных методов решения СЛАУ Алгоритм решения системы Ax = b итерационными методами включает следующие этапы. 1. Для применения итерационных методов система должна быть приведена к эквивалентному виду x = Bx + d. 2. Выбирается начальное приближение к решению системы уравнений x(0) = (x10, x20, …, xm0). 3. Находится последовательность приближений к корню. Для сходимости итерационного процесса достаточно, чтобы было выполнено условие . Критерий окончания итераций зависит от применяемого итерационного метода. 5.2. Метод Якоби Самый простой способ приведения системы к виду удобному для итерации состоит в следующем: из первого уравнения системы выразим неизвестное x1, из второго уравнения системы выразим x2, и т. д. В результате получим систему уравнений с матрицей B, в которой на главной диагонали стоят нулевые элементы, а остальные элементы вычисляются по формулам: , i, j = 1, 2, ... n. Компоненты вектора d вычисляются по формулам: , i = 1, 2, ... n. Расчетная формула метода простой итерации имеет вид , 35 или в покоординатной форме записи выглядит так: , i = 1, 2, ... m. Критерий окончания итераций в методе Якоби имеет вид: , где . Если , то можно применять более простой критерий окончания итераций: . Пример 1. Решение системы линейных уравнений методом Якоби Пусть дана система уравнений: Требуется найти решение системы с точностью  = 10-3 . Приведем систему к виду удобному для итерации: Выберем начальное приближение, например, - вектор правой части. Тогда первая итерация получается так: Аналогично получаются следующие приближения к решению. 36 , , Найдем норму матрицы . Будем использовать норму модулей элементов в каждой строке равна 0.2, то окончания итераций в этой задаче . Так как сумма = 0.2 < 1/2, поэтому критерий . Вычислим нормы разностей векторов: , Так как . , заданная точность достигнута на четвертой итерации. Ответ: x 1 = 1.102, x 2 = 0.991, x 3 = 1.101 5.3. Метод Зейделя Метод можно рассматривать как модификацию метода Якоби. Основная идея состоит в том, что при вычислении очередного (n+1)-го приближения к неизвестному xi при i >1 используют уже найденные (n+1)-е приближения к неизвестным x1, x2, ..., xi - 1, а не n-ое приближение, как в методе Якоби. Расчетная формула метода в покоординатной форме записи выглядит так: , i = 1, 2, ... m.. Условия сходимости и критерий окончания итераций можно взять такими же, как в методе Якоби. Пример 2. Решение систем линейных уравнений методом Зейделя Рассмотрим параллельно решение 3-х систем уравнений: , , Приведем системы к виду удобному для итераций: , Заметим, что условие сходимости , выполнено только для первой системы. 37 Вычислим 3 первых приближения к решению в каждом случае. 1-ая система: , , , . Итерационный процесс сходится. Точное решение здесь x 1 = 1.4, x 2 = 0.2. 2-ая система: , , , Итерационный процесс разошелся. Точное решение x 1 = 1, x 2 = 0.2. 3-я система: , , , Итерационный процесс зациклился. Точное решение x 1 = 1, Примечание. Для постройте чертеж. x 1 = 2. геометрической интерпретации полученных результатов Теорема о сходимости метода Пусть матрица системы уравнений A - симметричная и положительно определенная. Тогда при любом выборе начального приближения метод Зейделя сходится. Дополнительных условий на малость нормы некоторой матрицы здесь не накладывается. 5.4. Метод простой итерации Если A - симметричная и положительно определенная матрица, то систему уравнений часто приводят к эквивалентному виду: x = x - (Ax - b), где - итерационный параметр. Расчетная формула метода простой итерации в этом случае имеет вид: x (n+1) = x n - (Ax n - b). 38 Здесь B = E - A и параметр минимальной величину > 0 выбирают так, чтобы по возможности сделать . Пусть  min и  max - минимальное и максимальное собственные значения матрицы A. Оптимальным является выбор параметра . В этом случае принимает минимальное значение равное . ЗАДАЧИ ДЛЯ САМОСТОЯТЕЛЬНОГО РЕШЕНИЯ Задача 1. Убедиться в том, что если A- нижняя треугольная матрица с ненулевыми диагональными элементами, то метод Зейделя сходится за одну итерацию. Задача 2. Пусть A - верхняя треугольная матрица с ненулевыми диагональными элементами. Доказать, что метод Зейделя сходится за конечное число итераций и указать, за какое именно. Задача 3. Вывести оценку числа итераций, требуемых для достижения заданной точности в методе Якоби. КОНТРОЛЬНЫЕ ВОПРОСЫ Вопрос 1. Сформулируйте достаточное условие сходимости методов Якоби и метода Зейделя. Вопрос 2. Сформулируйте критерий окончания итераций в методе Якоби. Вопрос 3. Сформулируйте условия сходимости метода простой итерации и метода Зейделя для случая симметрических положительно определенных матриц. Вопрос 4. Из каких условий выбирается итерационный параметр в методе простой итерации. Вопрос 5. Сформулируйте алгоритм нахождения оптимального итерационного параметра в методе простой итерации. 39 Лекция 6 РЕШЕНИЕ СИСТЕМ НЕЛИНЕЙНЫХ УРАВНЕНИЙ Содержание      Системы нелинейных уравнений. Метод простой итерации. Метод Зейделя. Метод Ньютона-Рафсона. Методы контроля сходимости итерационных методов решения нелинейных систем 6.1. Системы нелинейных уравнений При решении задач моделирования на схемотехническом этапе проектирования ИС приходится решать системы уравнений, нелинейных по отношению к переменным. В отличие от систем линейных уравнений для систем нелинейных уравнений не известны прямые методы решения. Лишь в отдельных случаях систему можно решить непосредственно. Например, для системы из двух уравнений иногда удается выразить одно неизвестное через другое и таким образом свести задачу к решению одного нелинейного уравнения относительно одного неизвестного. Поэтому для нелинейных систем особую актуальность приобретают итерационные методы. Системы n линейных уравнений с n неизвестными x1, x2, ..., xn в общем случае принято записывать следующим образом: , где F1, F2 ,…, Fn – любые функции независимых переменных, в том числе и нелинейные (трансцендентные) относительно неизвестных. Как и в случае систем линейных уравнений, решением системы является такой вектор (или векторы) X*, который при подстановке обращает одновременно все уравнения системы в тождества. 40 . Система уравнений может не иметь решений, иметь единственное решение, конечное или бесконечное количество решений. Вопрос о количестве решений должен решаться для каждой конкретной задачи отдельно. Рассмотрим несколько простейших итерационных методов решения систем нелинейных уравнений, а именно, метод простой итерации, метод Зейделя и метод Ньютона. 6.2. Метод простой итерации Для реализации этого метода решаемую систему уравнений необходимо путем алгебраических преобразований привести к следующему виду, выразив из каждого уравнения по одной переменной следующим образом: . Выбирая затем вектор начального приближения , подставляют его в преобразованную систему уравнений. Из первого уравнения получают новое приближение к первой переменной, из второго – второй и т. д. Полученное уточненное значение переменных снова подставляют в эти уравнения и т.д. Таким образом, на (i+1)-м шаге итерационной процедуры имеем 41 . 6.3. Метод Зейделя Модификация Зейделя алгоритма простой итерации заключается в использовании уточненных значений переменных уже на текущем итерационном шаге. Так, для уточнения значений первой переменной используются только значения предыдущего шага, для второй переменной – значение x1 текущего шага, а остальных – от предыдущего и т.д.: 6.4. Метод Ньютона-Рафсона Математической основой метода является линеаризация функций F1, F2 ,…, Fn путем разложения в ряд Тейлора в окрестности точки начального приближения к решению и пренебрежением всеми членами ряда кроме линейных относительно приращений переменных. Для одной переменной ряд Тейлора в окрестности некоторой точки x = x0 выглядит так: . Для функции F1 от n переменных линейная часть разложения в ряд Тейлора в окрестности точки получается аналогичным образом: . 42 Введем следующие обозначения: xi = (xi - xi0) – приращение i-переменной – значение первой частной производной j-й функции по переменной xi при значении переменных (x10 , x20 ,…, xn0). Fj – значение j-й функции при значении переменных. (x10 , x20 ,…, xn0). После преобразования получим систему линейных уравнений порядка n относительно приращения переменных xi . . Или, в матричной форме, . В сокращенной матричной форме последнее выражение можно записать так: F' Δx = -F , где матрица значений частных производных F' носит название матрицы Якоби или якобиана системы уравнений. Решение этой системы дает вектор поправок к начальному приближению. Сложение его с вектором начального приближения дает новые, уточненные значения переменных: . Итерационная процедура далее продолжается аналогично. 43 6.4.1. Алгоритм решения системы нелинейных уравнений методом Ньютона -Рафсона 1. Система нелинейных уравнений приводится к стандартному виду и выбирается вектор начального приближения. 2. Рассчитывается матрица Якоби, членами которой являются значения частных производных в точке начального приближения 3. Решается система линейных уравнений относительно приращений (вектора поправок) переменных. 4. Находится очередное более точное значение вектора искомых переменных путем прибавления к вектору начального приближения вектора приращений 5. Проверяется условие сходимости к решению и, если оно не достигнуто, то процедура повторяется с пункта 2 . Частные производные, необходимые для расчета матрицы Якоби, можно рассчитать аналитически или же, если это невозможно или затруднительно, получать по формулам приближенного дифференцирования, например, как отношение приращения функции к приращению аргумента , где  – достаточно малое число. 6.5. Методы контроля сходимости итерационных методов решения нелинейных систем 1. Норма (эвклидова или максимум) вектора невязок или max (F1(i), F2(i) ,…, Fn(i))  . 2. Эвклидова норма вектора относительных отклонений переменных . 3. Норма-максимум вектора относительных отклонений 44 . КОНТРОЛЬНЫЕ ВОПРОСЫ Вопрос 1. Объясните сущность метода простой итерации. Вопрос 2. В чѐм отличие метода Зейделя от метода простой итерации?. Вопрос 3. Сформулируйте алгоритм решения систем нелинейных уравнений методом Ньютона-Рафсона. Вопрос 4. Сформулируйте критерии сходимости итерационных методов решения систем нелинейных уравнений. 45 Лекция 7 ОБРАБОТКА ДАННЫХ. ЧИСЛЕННЫЕ МЕТОДЫ ПРИБЛИЖЕНИЯ ФУНКЦИЙ Содержание          Области применения задачи приближения функций. Основные виды интерполяции, экстраполяция и аппроксимация. Постановка задач сглаживания функций. Задача интерполяции. Интерполяционный полином Лагранжа. Интерполяционный полином Ньютона. Погрешность полиномиальной интерполяции. Сплайн-интерполяция. Тригонометрическая интерполяция. 7.1. Области применения задачи приближения функций В технологии обработки данных важное место занимают методы приближения функций более простыми, хорошо изученными функциями, методы численного дифференцирования и численного интегрирования. При этом исследуемая приближаемая функция может быть задана как в аналитическом, так и дискретном виде (в виде экспериментальной таблицы). Подобная задача (иначе интерполяция функций) возникает в тех случаях, когда:  функция задана в виде таблицы, и необходимо знать значения функции для промежуточных значений аргументов, расположенных в таблице между узлами xi, j = 0,…. n, а также для аргументов, расположенных вне таблицы;  известна лишь таблица функции и требуется определить ее аналитическое выражение;  известно аналитическое выражение функции, но оно имеет очень сложный вид, вследствие чего возникает необходимость представления этой функции в более простом виде. Например, при вычислении определенных интегралов вида b  f  x dx a можно заменить подынтегральную функцию f(x) некоторой приближенной функцией P(x) в виде многочлена. Тогда b b a a  f x dx   Px dx . 46 7.2. Основные аппроксимация виды интерполяции, экстраполяция и Разберемся с рядом терминов, используемых в литературе при описании технологии приближения функций. Итак, пусть некоторая функция f(x) определена рядом своих узловых точек (xi, yi) на некотором отрезке [a, b]. Под интерполяцией мы будем подразумевать вычисление значений f(x) в любом промежутке [xi, xi+!] в пределах отрезка [a, b]. Соответственно, любое вычисление f(x) вне отрезка [a, b] является экстраполяцией. Наиболее распространены следующие виды интерполяции:  линейная интерполяция, при которой промежуточные точки, расположенные между двумя узловыми точками (xi, yi) и (xi+1, yi+1), лежат на отрезке прямой, соединяющей две ближайшие узловые точки;  квадратичная интерполяция, при которой промежуточные точки между узловыми точками (xi, yi), (xi+1, yi+1) и (xi+2, yi+2) лежат на отрезке параболы, соединяющей эти узловые точки;  полиномиальная интерполяция, при которой промежуточные точки вычисляются как значение некоторого многочлена pn(x), имеющего значения в узловых точках точно совпадающие с fi(xi);  Сплайновая интерполяция, при которой промежуточные точки находятся с помощью отрезков полиномов невысокой степени, проходящих через узловые точки и поддерживающие определенные условия стыковки в концевых точках; Подчеркнем ещѐ раз, что при интерполяции задача сводится к вычислению значений функции между узловыми точками (узлами), в то время как при экстраполяции — к вычислению функции вне того интервала, на котором она задана в виде таблицы, графически или иным способом. При аппроксимации таблично заданная функция (что, кстати, не является обязательным признаком аппроксимации) заменяется другой функцией, как правило, более простой и поэтому более быстро вычисляемой. Последняя приближенно описывает поведение исходной функции на некотором отрезке. При этом на различных отрезках аппроксимирующие функции могут быть (и чаще всего бывают) разными. 7.3. Постановка задачи сглаживания функций Пусть дана некоторая функция f(x) на отрезке x[a, b], которая является довольно сложной для исследования. Требуется заменить эту функцию некоторой простой, но хорошо исследуемой функцией (например, многочленом). Для этого на базе функции f(x) строят таблицу (ее называют сеточной функцией), которую можно заменить (сгладить) простой функцией с контролируемой погрешностью. Таблица 7.1. Сеточная функция, представленная в виде таблицы xi yi x0 y0 x1 y1 ... ... xn yn Рассмотрим два подхода к такой замене: 1. Пусть приближенная функция является многочленом n-й степени, 47 F(x) = Pn(x)=a0 + a1x+…+anxn (7.1) где n +1 - число узлов в табл. 7.1, с неизвестными параметрами ai, i = 0,…, n. Эта функция является приближением сеточной функции, т.е. yi = Pn(xi), i = 0,…, n. (7.2) В этом случае говорят, что функция (полином) Pn(x) интерполирует сеточную функцию, представленную табл. 7.1, а сама задача приближения называется задачей интерполяции. Точки xi, i = 0,…, n называют узлами интерполяции, а условие (7.2) – условием интерполяции. Введение интерполирующей функции позволяет вычислять значения f (x) не только в узлах интерполяции, но и между ними в точках ξ[xi-1, xi], i = 0,…, n, причем F(ξ) ≈ f(ξ). 2. При большом количестве точек xi, i = 0,…, n интерполяция требует большой гладкости (по n-й производной), что практически выполнить невозможно. Поэтому сглаживание сеточной функции (табл. 7.1) осуществляют путем минимизации некоторого функционала, построенного с помощью табл. 7.1 и многочлена 7.1 степени m, например, квадратичного функционала: (7.3) Процедуру сглаживания в этом случае называют аппроксимацией заданной функцией (7.1), в частности, аппроксимацию с использованием функционала (7.3) называют аппроксимацией с помощью точечного метода наименьших квадратов или просто метода наименьших квадратов (МНК). Если коэффициенты сглаживающей минимизации функционала функции (7.1) определяются путем (7.4) сглаживание называют интегральным методом наименьших квадратов. Если в качестве сглаживаемой функции задана экспериментальная таблица, то в методах сглаживания практически ничего не изменяется. Изменяются только методы оценки погрешности сглаживания. 7.4. Задача интерполяции Пусть на отрезке x[a, b] задана функция f(x), с помощью которой построена сеточная функция (табл. 7.1) или задана экспериментальная таблица (табл. 7.1). Определение. Точки x0 , …, xn называются узлами интерполяции. Требуется найти аналитическое выражение функции F(x), совпадающей в узлах интерполяции со значениями данной функции, т.е. F(x0) = y0, F(x1) = y1, …, F(xn) = yn . (7.5) Определение. Процесс вычисления значений функции F(x) в точках отличных от узлов интерполирования называется интерполированием функции f(x). Если 48 x[x0, xn], то задача вычисления приближенного значения функции в точке х называется интерполированием, иначе — экстраполированием. Геометрически задача интерполирования функции одной переменной означает построение кривой, проходящей через заданные точки (x0, y0), (x1, y1), …, (xn, yn) (рис. 7.1). В такой постановке задача может иметь бесконечное число решений или не иметь ни одного решения. Однако эта задача становится однозначно разрешимой, если вместо произвольной функции F(x) искать полином Pn ( x ) степени не выше n, удовлетворяющий условиям (7.5), т. е. такой, что При сглаживании функции (или экспериментальной таблицы) с помощью интерполяции в соответствии с условием интерполяции (7.2) значение интерполирующей функции и значение заданной функции в узлах сетки должны быть одинаковыми, следовательно, погрешность интерполяции в узлах xi, i = 0,…, n равна нулю (рис. 7.1). Рис. 7.1. Иллюстрация процедуры интерполяции Задача интерполяции имеет не единственное решение, но в случае, когда интерполирующей функцией является многочлен n-й степени (требующий n+1 узел интерполяции) вида (7.1) интерполяция имеет единственное решение, т.е. коэффициенты a0,…, an (в количестве n+1) определяются единственным образом. Действительно, используя табл. 7.1, можно составить систему из n+1 линейных уравнений (СЛАУ) относительно неизвестных коэффициентов a0,…, an: (7.5) 49 Матрица коэффициентов этой СЛАУ (7.5) носит специальное название — матрицы Вандермонда. Ее определитель не равен нулю, поскольку все значения узлов интерполяции различны между собой и ни одна из строк не является линейной комбинацией других строк, т.е. Таким образом, СЛАУ, а значит и задача полиномиальной интерполяции имеет единственное решение, поэтому вектор коэффициентов a0 ,…, an может быть выбран единственным образом. Примечание. Отметим, что вычисление коэффициентов полинома посредством решения системы (7.5) в вычислительной практике реализуется крайне редко. Причиной этого является плохая обусловленность матрицы Вандермонда, приводящая к заметному росту погрешности в выполнении условий интерполирования уже при сравнительно невысоких порядках полинома. К этому следует добавить, что вычислительные затраты реализации метода пропорциональны n3. 7.4.1. Интерполяционный полином Лагранжа Для полиномиальной интерполяции можно не решать СЛАУ (7.5), а записать многочлен (7.1) в специальной форме, одной из разновидностей которой является полином Лагранжа. Пусть функция f ( x) задана таблицей. Построим интерполяционный полином Ln(x), степени не больше n и для которого выполнены условия (7.6). Pn(x) = a0 + a1x + a2x2 + … + anxn , Pn(x0) = y0 , Pn(x1) = y1 , …, Pn(xn) = yn . (7.6) Будем искать Ln(x) в виде Ln ( x)  l0 ( x)  l1 ( x)  ...  ln ( x) (7.7) где li(x) – полином степени n , причем  y , если i  k li ( xk )   i  0, если i  k (7.8) Требование (7.8) совместно с (7.7) обеспечивает выполнение условий (7.6). Полиномы li(x) составим следующим образом: li ( x)  ci ( x  x0 )( x  x1 ) ... ( x  xi 1 )( x  xi 1 ) ... ( x  xn ) (7.9) Здесь ci – постоянный коэффициент, значение которого найдем из первой части условия (7.8): 50 ci  yi ( xi  x0 )( xi  x1 ) ... ( xi  xi 1 )( xi  xi 1 ) ... ( xi  xn ) Подставим ci в (7.9), тогда с учѐтом (7.7), интерполяционный полином Лагранжа будет иметь вид: n Ln ( x)   yi i 0 ( x  x0 )( x  x1 )...( x  xi 1 )( x  xi 1 )...( x  xn ) . ( xi  x0 )( xi  x1 )...( xi  xi 1 )( xi  xi 1 )...( xi  xn ) (7.10) Хотя формула Лагранжа имеет не совсем привычную форму, тем не менее, это одна из форм записи полинома n степени. Пример 7.1. Построить интерполяционный полином Лагранжа для функции, заданной таблицей: x f ( x) 1 12 3 4 4 6 Решение Из таблицы видно, что x0  1 , x1  3 , x2  4 , т. е. степень n интерполяционного полинома не выше, чем вторая. Используя формулу (7.10) получаем: L2 ( x)  12 ( x  3)( x  4) ( x  1)( x  4) ( x  1)( x  3) 4 6  2( x 2  7 x  12)  (1  3)(1  4) (3  1)(3  4) (4  1)(4  3) 2( x 2  5 x  4)  2( x 2  4 x  3)  2 x 2  12 x  22 Формуле Лагранжа (8) можно придать более компактный (сжатый) вид. Обозначим произведение: Пn1 ( x)  ( x  x0 )( x  x1 ) ... ( x  xn ) Продифференцируем это произведение Пn+1(x) по x: n Пn1 ( x)   ( x  x0 )...( x  xi 1 )( x  xi 1 )...( x  xn ) i 0 При x = xi (i = 0, 1, …, n) имеем: Пn1 ( xi )  ( xi  x0 ) ... ( xi  xi 1 )( xi  xi 1 ) ... ( xi  xn ) Тогда формула Лагранжа примет вид: n Ln ( x)   yi i 0 ( x  x0 )( x  xi 1 )( x  xi 1 )( x  xn )  ( xi  x0 )( xi  xi 1 )( xi  xi 1 )( xi  xn ) n Пn1 ( x) yi   yi  Пn1 ( x) ( x  xi ) Пn1 ( xi ) i 0 i 0 ( x  xi ) Пn 1 ( xi ) n Окончательно: 51 n yi ( x  xi )П n 1 ( xi ) Ln ( x )  П n 1 ( x ) i 0 (7.11) Полином (7.11) называют интерполяционным полиномом Лагранжа n-й степени, т.к. он, во-первых, удовлетворяет условию интерполяции Ln(xi) = yi , i = 0,…, n , и, во-вторых, имеет n-ю степень. С практической точки зрения интерес представляют два частных случая полинома Лагранжа:  При n = 1 имеем две точки, поэтому в данном случае формула Лагранжа представляет уравнение прямой y = L1(x), проходящей через две эти точки: y x b x a y0  y1 , ab ba где a и b — абсциссы этих точек.  При n = 2 получим уравнение параболы y = L2(x), проходящей через три точки ( x  b)( x  c ) ( x  a)( x  c ) ( x  a)( x  b) y y0  y1  y2 , (a  b)(a  c ) (b  a)(b  c ) (c  a)(c  b) где a, b и с — абсциссы этих точек. Формула Лагранжа удобна в задаче интерполирования многих функций в одной точке x, т. к. все значения множителей Lk(x) можно вычислить однажды для всех функций. Но формула Лагранжа имеет существенный недостаток. Бывает, что заданное число узлов недостаточно для достижения заданной точности. Тогда к данным узлам в табл. 7.1 добавляют еще один или несколько, что приводит к необходимости пересчета заново всех слагаемых в выражении 7.10 (или 7.11). Этого недостатка лишены интерполяционные полиномы Ньютона. 7.4.2. Интерполяционный полином Ньютона При построении интерполяционного полинома в форме Ньютона используется понятие разделенной разности, представляющее собой аналог понятия производной применительно к сеточным функциям. Разделенной разностью сеточной функции (табл. 7.1) нулевого порядка в узлах xi , i = 0,…, n называются значения этой функции в этих узлах f (xi)= yi , i = 0,…, n. Определение 1. Разделенной разностью сеточной функции первого порядка в узлах xi , i = 0,…, n-1 называют отношение Определение 2. Разделенной разностью сеточной функции второго порядка в узлах xi , i = 0,…, n-2 называют отношение 52 Определение 3. Разделенной разностью сеточной функции n-го порядка в узле x0 называют отношение С использованием разделенных разностей интерполяционный полином Ньютона записывается в форме: (7.12) Отметим, что при добавлении новых узлов первые члены многочлена Ньютона остаются неизменными. Если функция задана в точках x0, x1,…, xn, то при построении интерполяционного многочлена Ньютона удобно пользоваться таблицей, называемой таблицей разделенных разностей, пример которой для n = 4 приведен в табл. 7.2. Таблица 7.2 Для повышения точности интерполяции в сумму (7.7) могут быть добавлены новые члены, что требует подключения дополнительных интерполяционных узлов. При этом безразлично, в каком порядке подключаются новые узлы. Этим формула Ньютона выгодно отличается от формулы Лагранжа. 7.4.3. Погрешность полиномиальной интерполяции Ясно, что в узлах интерполяции погрешность интерполяционного полинома Ln(x) или Nn (x) равна нулю Погрешность Ln(x) - f(x) представляющая собой разность между значением интерполяционного многочлена Ln(x) и значением функции f(x) в точке x, не совпадающей с узлом интерполяции имеет вид: 53 (7.13) Где - точка, в которой ищется погрешность (не совпадает с узлами интерполяции). Поскольку точка ξ(a,b) неизвестна, то вместо погрешности (7.13) вводится верхняя оценка погрешности в виде (7.14) которая и используется на практике. Таким образом, погрешность интерполяции зависит как от величины соответствующей производной приближаемой функции, так и от расположения узлов. Минимизировать погрешность приближения достаточно гладкой функции на отрезке [a,b] полиномом можно, расположив узлы интерполяции следующим образом: где - корни полинома Чебышева Hn (x) = cos(n arccos(x)) или в рекуррентном виде Отметим также, что такое расположение узлов интерполяции гарантирует сходимость интерполяционного полинома к приближаемой функции при повышении числа узла интерполяции (степени полинома), тогда как при равномерном распределении узлов в ряде случаев может наблюдаться расходимость (такая ситуация хорошо иллюстрируется известным примером Рунге, в котором функция приближается интерполяционным полиномом на отрезке [-1,1]). Пример 7.2. Построить интерполяционный полином Лагранжа, совпадающий с функцией f(x) = 3x, x[-1,1] в точках x0 = -1, x1 = 0, x2 = 1. Вычислить значение сеточной функции и оценить погрешность интерполяции в точке x* = 5,0. Решение Составим сеточную функцию и занесем ее в таблицу. Поскольку n=2, то необходимо построить интерполяционный полином L2(x). 54 Проверим условия интерполяции Значение сеточной функции в точке x* = 0,5 вычислим по интерполяционному многочлену y(0,5) = L2(0,5) = 1,8333. Верхнюю оценку погрешности интерполяционного многочлена определим в соответствии с выражением (7.14) Поскольку функция f(x) = 3x известна, то можно вычислить точное значение абсолютной погрешности в точке x*= 0,5. Пример 7.3. Составить интерполяционный полином Ньютона для заданной таблицы Решение Таблица задана с неравномерным шагом, поэтому для решения задачи воспользуемся многочленом Ньютона с разделенными разностями (формула (7.12)) где f(x0) = y0 = 1/3; 55 Таким образом, Условия интерполяции соблюдены: 7.4.4. Сплайн-интерполяция Рассмотренная в предыдущих разделах интерполяция, когда интерполяционный полином строится сразу по всем узлам интерполяции называется глобальной интерполяцией. При этом увеличение числа узлов автоматически приводит к повышению степени полинома, и как следствие, к проявлению его колебательных свойств (рис. 7.2). Поэтому обычную полиномиальную интерполяцию осуществляют максимум по 3-4 узлам. Интерполяцию по нескольким узлам табл. 7.1, называют локальной. Рис. 7.2. Проявление колебательных свойств глобального интерполяционного полинома Рис. 7.3. Локальная интерполяция по каждым трем узлам Однако, такая локальная интерполяция с помощью Ln или Nn страдает тем недостатком, что интерполирующая функция в узлах стыковки полиномов имеет 56 непрерывность только нулевого порядка, т.е. интерполирующая принадлежит классу функций C0 (см. рис. 7.3 для Ln и Nn в узле x*). функция От этих недостатков свободна сплайн-интерполяция, которая требует непрерывности в узлах стыковки локальных многочленов по производным соответственно порядка один, два, и т.д. Определение. Сплайном - степени m дефекта r называется (m - r) раз непрерывно дифференцируемая функция, которая на каждом отрезке [xi – xi-1], i = 1,…, n представляет собой многочлен степени m. Наиболее распространенными в науке и технике являются сплайны 3-й степени дефекта один, т.е. т.е. дважды непрерывно дифференцируемый многочлен 3-й степени на каждом отрезке [xi , xi-1], i = 1,…, n. Сплайны, удовлетворяющие условию интерполяции, называются интерполяционными. Основным достоинством интерполяционного кубического сплайна дефекта один является следующее: этот сплайн обладает минимумом интегральной кривизны на всем заданном отрезке [a, b] по сравнению с другими интерполяционными функциями f *(x), т.е. . Геометрически это означает, что если тяжелую упругую нить повесить на ряд гвоздей, то она примет форму кубического сплайна дефекта 1, приведенную на рис. 7.4. Рис. 7.4. Тяжелая упругая нить, геометрически представляющая собой кубические сплайны дефекта один Рассмотрим алгоритм построения интерполяционного кубического сплайна S(x), i =1,…, n дефекта один в соответствии с табл. 7.2. Кубический полином Si(x) на отрезке [xi-1 , xi] имеет четыре неизвестных коэффициента. Количество отрезков в соответствии с табл. 7.2 равно n. Для определения 4×n коэффициентов можно использовать следующие условия в узлах интерполяции:  условие интерполяции S(xi) = yi для i =1,…, n  непрерывность сплайна S(xi - 0) = S(xi + 0) для i =1,…, n-1 57  непрерывность производных 1-го порядка S’(xi - 0) = S’(xi + 0) для i =1,…, n-1  непрерывность производных 2-го порядка S”(xi - 0) = S”(xi + 0) для i =1,…, n-1 Таким образом, всего имеется (n+1) + 3(n-1) = 4n - 2 условий. В качестве двух недостающих условий задают значения производных 1-го или 2-го порядка в узлах x0 и xn. Для вывода используем значения S”(x0) = S”(xn) . В этом случае сплайн называется естественным. Пусть S”(x) = q(x). Рассмотрим поведение функции q(x) на отрезке [xi-1 , xi] (см. рис. 7.5). Рис. 7.5. Поведение функций S”(x) на элементарных отрезках Поскольку сплайн является многочленом 3-й степени, то на каждом отрезке [xi-1 , xi], 2-я производная будет линейна. Найдем ее с помощью интерполяционного многочлена Лагранжа 1-й степени L1(x): (7.15) Выражение (7.15) уже удовлетворяет условиям непрерывности производных 2-го порядка. Действительно, подставим в (7.15) x = xi − 0, получим q(xi − 0) = qi. Затем, записывая выражение (7.15) для отрезка x[xi , xi+1] и подставляя в него xi + 0 вместо x, получим q(xi + 0) = qi, что и требовалось показать. Для нахождения сплайна проинтегрируем дважды выражение (7.10), получим (7.16) где C1 и C2 найдем из удовлетворения значений сплайна (7.16) в узлах xi-1, xi условиям интерполяции 58 . Решая эту СЛАУ относительно C1 и C2 и подставляя их в (7.16), найдем следующее выражение для сплайна степени 3 дефекта 1: (7.17) В этом сплайне узловые значения для вторых производных qi пока неизвестны. Будем искать их из условий непрерывности первых производных в узлах xi. Для нахождения производной S′(xi + 0) запишем (7.17) для отрезка [xi , xi+1] (7.18) Вычисляя от (7.17) и (7.18) производные первого порядка и подставляя в них значение x = xi, получим Приравняем эти выражения в соответствии с условиями непрерывности первых производных в узлах интерполяции xi, получим (7.19) q0 = qn = 0. (7.20) Система (7.19) с заданными краевыми условиями (7.20) представляет собой СЛАУ относительно qi = S’’(xi), i = 1,…, n -1, имеет трехдиагональную матрицу и, следовательно, ее можно решать методом прогонки. Подставляя найденные qi, i = 0,…, n, в (7.18), получим кубические сплайны дефекта один на каждом отрезке x[xi , xi+1], i = 1,…, n. Таким образом, определяющими выражениями для определения кубических сплайнов дефекта один являются выражения (7.18), (7.19), (7.20). Пример 7.4. Для заданной таблицы с шагом h = xi - xi-1 = 1 = const построить интерполяционный кубический сплайн дефекта один, выписав соответствующие 59 уравнения на каждом отрезке x[xi-1 , xi], i = 1,…, 4. Проверить непрерывность сплайна и его производных до второго порядка включительно в узле x* = 2. Решение Под заданной таблицей сформируем дополнительную строку для вторых производных сплайнов S’’(xi) ≡ qi , которая заполняется по мере их вычисления (сразу можно вписать в нее q0 = q4 = 0). Для узлов x1= 2; x2 = 3; x3 = 4 с учетом q0 = q4 = 0 составляется СЛАУ (7.19) относительно неизвестных q1, q2, q3. Вычисляются прогоночные коэффициенты по формулам и значения Заносим эти значения в дополнительную строку таблицы и для каждого из четырех интервалов выписываем уравнения сплайна (7.12). 60 Проверим правильность построения сплайна для узла x* = 2. К данному узлу примыкают кривые S1(x) и S2(x). 7.2.5. Тригонометрическая интерполяция Поскольку многие явления в природе имеют циклический характер, широкое практическое применение получила интерполяция дискретных периодических функций тригонометрическими полиномами вида: где - частота k-ой гармоники, L - период, ak, bk - коэффициенты разложения. Такой подход позволяет представить сложную циклическую структуру в виде суперпозиции простых периодических функций (элементарных гармоник). 61 Рассмотрим таблично заданную на периоде L функцию yi(xi) с равномерным распределением узлов xi = x0 + ih , h=L/n , i = 0,…, n. Тогда, если n - четно (n = 2m), существует единственный интерполяционный тригонометрический полином Tm(x) степени m = n /2, удовлетворяющий условиям Tm(xi) = yi, i = 0,…, n: Коэффициенты разложения определяются следующим образом: Отметим, что периодичность исходной функции yi(xi) предполагает, что y0 = yn. Если это условие не выполнено, то построенный тригонометрический полином будет удовлетворять условиям интерполяции во всех узлах, кроме последнего, т.е. Tm(xi) = yi, i = 0,…, n – 1. В последнем узле будет выполняться условие периодичности Tm(xn) = Tm(x0). КОНТРОЛЬНЫЕ ВОПРОСЫ Вопрос 1. Когда возникает необходимость в использовании интерполяционных методов? Вопрос 2. В чѐм сущность задачи интерполирования? Вопрос 3. Поясните смысл терминов: интерполяция, экстраполяция и аппроксимация. Вопрос 4. Как строится интерполяционный многочлен Лагранжа? Вопрос 5. Как оценивается погрешность интерполяционной формулы Лагранжа? Вопрос 6. Дайте определение понятий разделенной разности нулевого и первого порядков. Вопрос 7. Объясните принцип построения интерполяционного полинома Ньютона. Вопрос 8. В чем разница между глобальной и локальной разновидностями интерполяции? Вопрос 9. Что такое сплайн-интерполяция? 62 Лекция 8 ЧИСЛЕННЫЕ МЕТОДЫ ВЫЧИСЛЕНИЯ ОПРЕДЕЛЕННЫХ ИНТЕГРАЛОВ Содержание  Сущность численных методов интегрирования.  Использование квадратурных формул  Интегрирование методом прямоугольников  Интегрирование методом трапеций  Метод Симпсона  Метод Монте-Карло 8.1. Сущность численных методов интегрирования Если функция f(x) непрерывна на отрезке [a, b] и известна ее первообразная F(x), то определенный интеграл от этой функции в пределах от a до b может быть вычислен по формуле Ньютона — Лейбница: b  f ( x)dx  F(b)  F(a) , a где F (x) = f(x). Однако во многих случаях функция F(x) не может быть найдена с помощью элементарных средств или является слишком сложной, вследствие этого вычисление определенного интеграла по приведенной формуле может быть практически невыполнимо. Кроме того, на практике подынтегральная функция f(x) часто задается таблично и тогда само понятие первообразной теряет смысл. Поэтому особое значение приобретают численные (приближенные) методы вычисления определенных интегралов. Особенно это актуально при решении задач научного и инженерно-технического характера. Общим способом интегрирования любых функций является численное интегрирование, методы которого в большинстве своем просты и легко переводятся на алгоритмические языки. Геометрически интеграл функции f(x) в пределах от a до b представляет собой площадь криволинейной трапеции (рис. 8.1), ограниченной графиком этой функции, осью x и прямыми x = a и x = b. Численные методы интегрирования используют замену площади криволинейной трапеции на конечную сумму площадей более простых геометрических фигур, которые могут быть вычислены точно. В этом смысле говорят об использовании квадратурных формул (по аналогии с задачей о квадратуре круга – построение квадрата с площадью, равной площади круга с определенным радиусом). 63 Рис.8.1. Графическое представление определенного интеграла от функции f(x) в пределах от x = a до b. Взятие интеграла эквивалентно вычислению площади под кривой 8.2. Использование квадратурных формул В большинстве методов используется приближенной представление интеграла в виде конечной суммы (квадратурная формула): , где ci – постоянные, называемые весами, xi – значения аргумента, принадлежащие интервалу [a, b] — называемые узлами. В основе квадратурных формул лежит идея аппроксимации на отрезке интегрирования графика подынтегрального выражения функциями более простого вида, которые легко могут быть проинтегрированы аналитически и, таким образом, легко вычислены. Наиболее просто задача построения квадратурных формул реализуется для полиномиальных математических моделей, рассмотренных ранее в лекции 7 «Обработка данных. Численные методы приближения функций». Напомним некоторые сведения из предыдущей лекции. Многочлен (полином) порядка n имеет вид , и определяется, таким образом, (n+1) значениями констант ai . Если известно значение функции в (n+1) точках i = 0, 1, …, n, то данные параметры легко определяются из системы (n+1) линейных уравнений с (n+1) переменными ai 64 . Если все xi различны, то данная система уравнений имеет единственное решение, так как определитель системы, составленный из коэффициентов системы линейных уравнений (так называемый определитель Вандермонда) будет отличен от нуля . Определив коэффициенты интерполяционного многочлена, можно легко вычислить приближенное значение интеграла, заменив подынтегральную функцию на полученный многочлен . . Узлы интерполирования на отрезке интегрирования могут быть расположены на равном удалении друг от друга (эквидистантные узлы). В этом случае расстояние между узлами для полинома степени n определится с помощью следующей формулы: , где h – шаг, xi – узлы интерполирования. В случае n = 0 получаем метод прямоугольников. . Здесь график функции f(x) на отрезке интегрирования [a, b] заменяется на горизонтальную линию (полином степени 0). 8.2.1. Интегрирование методом прямоугольников (метод Эйлера) Пусть функцию (рис. 8.1) необходимо проинтегрировать численным методом на отрезке [a, b]. Разделим отрезок на N равных интервалов (на рисунке N = 4). 65 Рис. 8.2. Схема численного расчета интеграла методом прямоугольника Площадь каждой из 4-х криволинейных трапеций можно заменить на площадь соответствующего прямоугольника. Ширина всех прямоугольников одинакова и равна . В качестве высоты прямоугольников можно взять значение функции на левой границе. В этом случае высота первого прямоугольника составит f(a), второго – f(x1), третьего – f(x2), последнего – f(x3). Приближенное значение интеграла получается суммированием площадей прямоугольников . Если в качестве выбора высоты прямоугольников взять значение функции на правой границе, то в этом случае высота первого прямоугольника составит f(x1), второго – f(x2), третьего – f(x3), последнего – f(b). . Очевидно, что первый способ дает приближение к интегралу с недостатком, а второй — c избытком. Можно предложить еще один способ более точный, чем обе эти формулы – использовать для аппроксимации значение функции в середине отрезка интегрирования. . В общем виде, если отрезок [a, b] разбить на N равных интервалов интегрирования и к каждому интервалу применить формулу прямоугольников, то для каждого из трех рассмотренных способов получим следующие обобщенные формулы: 66 , и . Погрешность метода. Всякое приближенное вычисление имеет определенную ценность лишь тогда, когда оно сопровождается оценкой допущенной при этом погрешности. Поэтому формулы прямоугольников будут практически пригодны для приближенного вычисления интегралов лишь в том случае, если будет существовать удобный способ оценки получающейся при этом погрешности (при заданном n), позволяющий к тому же находить и число частей n разбиения сегмента, гарантирующее требуемую степень точности приближенного вычисления. Будем предполагать, что функция f(x) имеет ограниченную производную на сегменте [a, b], так что существует такое число М>0, что для всех значений х из [a, b] выполняется неравенство |f'(x)|  M. Качественный смысл этого неравенства заключается в том, что скорость изменения значения функции ограничена. В реальных физических системах это требование практически всегда выполнено! Тогда погрешность Rn вычисления интеграла по формулам прямоугольников может быть оценена по формуле: |Rn|  M(b-a)2/2n При неограниченном возрастании n выражение M(b-a)2/2n, а следовательно, и погрешность Rn стремится к нулю, т.е. точность приближения будет тем больше, чем на большее число частей будет разделен сегмент [a, b]. Следовательно, для вычисления интеграла с заданной точностью  достаточно сегмент [a, b] разбить на число частей n = M(b-a)2/2 . Метод прямоугольников – это наиболее простой и вместе с тем наиболее грубый метод приближенного интегрирования. 8.2.2. Интегрирование методом трапеций Использование для интерполяции полинома первой степени (прямая линия, проведенная через две точки) приводит к формуле трапеций (рис. 8.3). Рис. 8.3. Принцип численного интегрирования методом трапеций 67 В качестве узлов интерполирования берутся концы отрезка интегрирования. Таким образом, криволинейная трапеция заменяется на обычную трапецию, площадь которой может быть найдена как произведение полусуммы оснований на высоту. . В случае N отрезков интегрирования площадь криволинейной трапеции аппроксимируется суммой из N обычных трапеций (рис. 8.4). Рис. 8.4. Графическая интерпретация численного интегрирования методом трапеций Здесь для всех узлов, за исключением крайних точек отрезка, значение функции войдет в общую сумму дважды (так как соседние трапеции имеют одну общую сторону). . Интересно, что формула трапеций может быть получена, если взять половину суммы формул прямоугольников по правому и левому краям отрезка . 68 Проиллюстрируем использование формулы трапеций на примере, представленном на рис. 8.2. . (8.1) Погрешность метода. Для оценки погрешности Rn, возникающей при замене исходной функции аппроксимирующей ее функцией (8.1), доказывается, что абсолютная величина погрешности удовлетворяет неравенству: , где М2 – максимум модуля второй производной функции на отрезке [a,b], т.е. Погрешность Rn будет меньше заданного числа   0, если взять . Таким образом, в методе трапеций ошибка вычисления интеграла приблизительно пропорциональна квадрату шага интегрирования h 2, т.е. Rn Rn ~ h 2 . Резюме Для вычисления интеграла некоторой функции в пределах a, b методом прямоугольника либо трапеций необходимо разделить отрезок [a, b] на nо интервалов и найти сумму их площадей (отдельных прмоугольников либо трапеций). Затем нужно увеличить число интервалов (n1), опять вычислить сумму элементарных площадей и сравнить полученное значение с предыдущим результатом. Это следует повторять до тех пор (ni), пока не будет достигнута заданная точность результата (критерия сходимости). Главное преимущество метода трапеций – простота. Однако если при вычислении интеграла требуется высокая точность, применение этого метода может потребовать слишком большого количества итераций или машинного времени. Пример 8.1. Вычислить интеграл численным методом, пользуясь правилом трапеций Решение Для n = 1 . 69 Для n = 2 . Для n = 4 . Для n = 64 h = 0,0156 I = 0,333 (точное решение равно 1/3). 8.2.3. Метод Симпсона Использование трех точек для интерполирования подынтегрального выражения позволяет использовать параболическую функцию (полином второй степени). Это приводит к формуле Симпсона приближенного вычисления интеграла (рис. 8.5, b). Рис. 8.5. Сопоставление принципов интегрирования методом трапеций (a) и параболой (b) Рассмотрим произвольный интеграл . Воспользуемся заменой переменной таким образом, чтобы границы отрезка интегрирования вместо [a,b] стали [-1,1], для этого введем переменную z: . Тогда 70 , и . Рассмотрим задачу интерполирования полиномом второй степени (параболой) подынтегральной функции, используя в качестве узлов три равноудаленные узловые точки: z = -1, z = 0, z = +1 (шаг равен 1, длина отрезка интегрирования равна 2). Обозначим соответствующие значения подынтегральной функции в узлах интерполяции . Система уравнений для нахождения коэффициентов полинома , проходящего через три точки: (-1, f-1), (0, f0) и (+1, f+1), примет вид или . Решая последнюю систему уравнений относительно коэффициентов. получаем . Вычислим теперь значение интеграла от интерполяционного многочлена . Путем обратной замены переменной вернемся к исходному интегралу. Учтем, что z = -1, соответствует , 71 z = 0, соответствует , z = +1, соответствует . В результате получим формулу Симпсона для произвольного интервала интегрирования: , где . При необходимости, исходный отрезок интегрирования может быть разбит на N сдвоенных отрезков, к каждому из которых применяется формула Симпсона. Шаг интерполирования при этом составит . Для первого отрезка интегрирования узлами интерполирования будут являться точки a, a+h, a+2h, для второго – a+2h, a+3h, a+4h, третьего a+4h, a+5h, a+6h и т.д. Приближенное значение интеграла получается суммированием N площадей: . В данную сумму входят одинаковые слагаемые (для внутренних узлов с четным значением индекса - 2i). Поэтому можно перегруппировать слагаемые в этой сумме таким образом , что эквивалентно , так как . Погрешность метода. Погрешность этого приближенного метода уменьшается пропорционально длине шага интегрирования в четвертой степени, т.е. при увеличении числа интервалов вдвое ошибка уменьшается в 16 раз 72 Rn ~ h 4 . Пример 8.2. Вычислить интеграл численным методом, пользуясь правилом Симпсона . Решение Для n = 1 . Для n = 2 . Примечание. Точное решение равно 0,2 8.2.4. Об оценке точности квадратурных формул Точность квадратурной формулы при фиксированном числе узлов существенно зависит от расположения этих узлов. При неудачном расположении узлов квадратурная формула может дать сильно искаженные результаты. При наличии значительного числа нулей подынтегральной функции или при большом числе ее экстремумов (т. е. когда много нулей производной f΄(x)) точность квадратурных формул сильно понижается. Поэтому шаг h следует выбирать так, чтобы он был намного меньше расстояний между соседними нулями (корнями) функции f(x) и ее производной f΄(x). Для этого рекомендуется, например, основной отрезок интегрирования [a,b] разбить на частичные отрезки [α,β], внутри каждого из которых, функции f(x) и f΄(x) сохраняют постоянный знак (если это возможно!), и производить вычисление интеграла по частям, выбирая для каждого отрезка, вообще говоря, свой шаг. 8.3. Метод Монте- Карло В некоторых случаях из-за особенности подынтегральной функции (например, изза ее большой сложности, неявном способе задания и т.д.), описанные выше методы нельзя или нецелесообразно использовать. В задачах, не требующих высокой точности, широкое распространение получил метод Монте-Карло. 73 Проиллюстрируем применение этого вычисления следующего интеграла: метода на примере приближенного . График подынтегрального выражения приведен на рис. 8.6. Очевидно, что точное значение интеграла равно четверти площади круга единичного радиуса. Рис. 8.6. Иллюстрация метода Монте-Карло Для знакомства с этим методом построим прямоугольную область, которая будет полностью включать в себя искомый интеграл. В данном случае это будет квадрат с единичным ребром, показанный на рис. 8.6. Далее, с помощью датчика случайных чисел генерируются точки (xi , yi) для i = 1, 2, 3,…, N, попадающие в эту область. В данном случае абсциссы и ординаты точек должны быть случайными числами, равномерно распределенными на отрезке [0, 1]. Для каждой точки проверяется, попадает ли она в область под или над графиком функции, то есть проверяется условие: . Если условие выполняется, то выбранная точка соответствует «успеху», если нет – то «промаху». Таким образом, процедура может быть описана как игра в «попадание» случайно выбранной точки в область под графиком (отсюда и название метода — МонтеКарло). Вполне очевидно, что отношение числа «попаданий» (Nусп) к общему числу попыток (N) должно в пределе стремиться к доли площади прямоугольной области (Sпр), которую занимает область под интегрируемой функцией (значение интеграла, I). . Отсюда получается формула метода Монте-Карло: . 74 Для реализации метода существенное значение имеет качество используемого датчика случайных чисел. Идеальный датчик должен давать равномерное распределение чисел в заданном диапазоне. Точность расчета интеграла определяется так же числом точек (испытаний) N, используемых при вычислениях и, очевидно, должна увеличиваться при его росте. Примечание. В лабораторной работе 8 «Численное интегрирование» используется модифицированный вариант рассмотренного здесь алгоритм МонтеКарло. Резюме Метод Монте-Карло требует значительных временных ресурсов и характеризуется медленной сходимостью к истинному решению. Однако в сложных ситуациях, когда функцию невозможно рассчитать аналитическими методами или при расчете многомерных интегралов, применение этого метода может стать единственно возможным вариантом. КОНТРОЛЬНЫЕ ВОПРОСЫ Вопрос 1. В каких случаях используются приближенные методы вычисления определенных интегралов? Вопрос 2. В чем сущность численного интегрирования? Вопрос 3. Какие квадратурные формулы вы знаете? Вопрос 4. Как выбирается шаг интегрирования? Вопрос 5. На чем основан метод Монте-Карло? 75 Лекция 9 ЧИСЛЕННОЕ ДИФФЕРЕНЦИРОВАНИЕ Содержание  Сущность численных методов дифференцирования.  Дифференцирование путем интерполяции полиномом  Дифференцирование методом конечных разностей  Оценка погрешности методом Рунге 9.1. Сущность численных методов дифференцирования К численному (приближенному) дифференцированию чаще всего прибегают, когда приходится вычислять производные от функций, заданных таблично, или, когда непосредственное дифференцирование затруднительно. Пусть в некоторой точке требуется вычислить производные первого, второго и т.д. порядков от дискретно заданной (сеточной) функции (табл. 9.1): Таблица 9.1. Таблично заданная функция xi yi x0 y0 x1 y1 ... ... xn yn При этом возможны два случая: а) точка x* (xi-1, xi) i = 1,…, n; б) точка x*= xi , i = 0,…, n, т.е. совпадает с одним из узлов заданной таблицы. 9.2. Интерполяция полиномом В первом случае заданная таблица сглаживается какой-либо функцией P(x), являющейся глобальным (локальным) интерполяционным полиномом, или полиномом, полученным с использованием МНК (метода наименьших квадратов) с некоторой погрешностью Rn(x), в результате чего имеют место следующие равенства: f(x) = P(x) + Rn(x), f(x*) = P(x*) + Rn(x*): f′(x) = P′(x) + R′n(x), f′(x*) = P′(x*) + R′n(x*): f′′(x) = P′′(x) + R′′n(x), f′′(x*) = P′′(x*) + R′′n(x*): …………………………………………………………. (9.1) Отметим, что численное дифференцирование представляет собой операцию менее точную (иногда говорят — некорректную), чем интерполирование. Действительно, близость друг другу ординат двух кривых y=f(x) и y=P(x) на отрезке [a,b] еще не гарантирует близости на этом отрезке их производных f′(x) и P′(x), т. е. малого расхождения угловых коэффициентов касательных к рассматриваемым кривым при одинаковых значениях аргумента (рис. 9.1). 76 Рис. 9.1. К понятию некорректности численного дифференцирования 9.3. Интерполяция конечными разностями Во втором случае (x*= xi , i = 0,…, n) используется аппарат разложения функций в ряд Тейлора, для чего функция в точке x* должна иметь достаточное число производных. С этой целью предполагается, что заданная таблица (9.1) является сеточной функцией для некоторой функции y(x) (т.е. yi = y(xi )), имеющей в точке производные до четвертого порядка включительно. Пусть внутренний узел x*= xi , i = 1,…, n-1 окружают другие узлы (рис. 9.2), причем xi+1 - xi = xi - xi-1 = h. Рис. 9.2. Схема численного дифференцирования методом конечных разностей В этом случае, разлагая значения yi-1 , yi+1 на точной функции y(x) в ряд Тейлора в окрестности точки xi до производной четвертого порядка включительно, получим: , (9.2) (9.3) Выразим вначале yi’ из (9.2), а затем из (9.3), разделив предварительно на h и оставляя слагаемые с первой степенью шага h, получим; 77 , (9.4) . (9.5) Вычтем из (9.3) выражение (9.2), разделим полученное соотношение на 2h, получим следующее значение производной первого порядка в точке xi (слагаемые с производными четного порядка сокращаются за исключением слагаемого с производной четвертого порядка): , (9.6) где — центральная разность первого порядка. Определим порядок аппроксимации (порядок точности) формулы численного дифференцирования производной как показатель степени h в главном члене погрешности. Выражение (9.4) определяет производную первого порядка с помощью отношения конечных разностей слева. При этом имеет место первый порядок аппроксимации относительно шага h. Выражение (9.5) определяет производную первого порядка с помощью отношения конечных разностей справа. При этом имеет место первый порядок аппроксимации относительно шага h. Выражение (9.6) определяет производную первого порядка с помощью отношения центральных разностей. В этом случае реализуется второй порядок аппроксимации относительно шага h. Таким образом, односторонние производные (9.4) и (9.5), имеющие 1-й порядок аппроксимации менее точны, чем центрально-разностная производная (9.6), имеющая второй порядок. Сложим выражения (9.2) и (9.3) (слагаемые с производными нечетного порядка сокращаются) и разделим на h2, получим: . (9.7) Выражение (9.7) определяет производную второго порядка в узле с помощью отношения центральных разностей второго порядка. Она имеет второй порядок аппроксимации относительно шага h. 9.4. Метод Рунге оценки погрешности и уточнения формул численного дифференцирования Из формул (9.4)-(9.7) видно, что метод p-го порядка численного дифференцирования совпадает с показателем степени шага h в главном члене погрешности и имеет вид: 78 (9.8) где а остаточный член имеет вид С целью оценки погрешности продифференцируем численно методом p-го порядка функцию f(xi) = yi , i = 0,…, n с шагом h, получим выражение (9.8). Затем продифференцируем численно функцию тем же методом p-го порядка, с шагом kh (k=1/2; 1/4; 1/16; ...), получим: . (9.9) Вычитая из (9.9) выражение (9.8) и определяя из полученного равенства ψ(x), находим: . (9.10) Выражение (9.10) можно использовать для оценки погрешности численного дифференцирования. Подставляя (9.10) в (9.9), получаем окончательно: . (9.11) Из последнего выражения видно, что это уже метод порядка p+1, т.е. на порядок точнее. КОНТРОЛЬНЫЕ ВОПРОСЫ Вопрос 1. В каких случаях используются приближенные методы вычисления производных? Вопрос 2. В чем сущность численного дифференцирования? Вопрос 3. Плюсы и минусы численного дифференцирования с помощью интерполяции исходной функции полиномом. Вопрос 4. Поясните сущность численного дифференцирования методом конечных разностей? Вопрос 5. Методика Рунге оценки погрешности. 79 Лекция 10 МЕТОДЫ РЕШЕНИЯ ОБЫКНОВЕННЫХ ДИФФЕРЕНЦИАЛЬНЫХ УРАВНЕНИЙ Содержание  Задача Коши.  Сущность процедуры численного решения ОДУ.  Одношаговые конечно-разностные методы решения задачи Коши.  Метод Эйлера.  Модифицированный метод Эйлера.  Метод Рунге-Кутта.  Метод Адамса. 10.1. Введение Уравнения, содержащие дифференциальными. одну или несколько производных, называются Инженерные и научные задачи часто связаны с решением дифференциальных уравнений, так как с помощью последних описываются многие физические явления. Соответственно процессы в технических устройствах так же описываются дифференциальными уравнениями. Природа этих процессов различна. При анализе тепловых режимов аппаратуры рассчитывают тепловые потоки, при изучении электромагнитных процессов – электрические и магнитные поля, при оценке прочности изделий вычисляют механические напряжения и деформации. Особую значимость эта тема играет при расчете переходных процессов электронных и интегральных схем. Практически все современные системы моделирования и САПР, включая различные клоны SPICE, базируются на мощных программных средствах численного решения систем дифференциальных уравнений. Более того, именно потребность в создании этих систем, послужило мощным стимулом для разработки и совершенствования данного класса численных методов. К сожалению, для многих практически важных случаев задачи, описываемые дифференциальными уравнениями, весьма сложны, и получить их точное решение оказывается затруднительно или невозможно. Эти трудности могут быть связаны с видом уравнения, например, с его нелинейным характером. Однако решить подобные сложные задачи также как и более простые можно с помощью компьютера. Поэтому методы решения дифференциальных уравнений на ЭВМ широко применяются в инженерной практике. Методы решения задач, содержащих обыкновенные дифференциальные уравнения, зависят от их математической формулировки. Рассмотрим их. 80 10.2. Классификация уравнений Дифференциальные уравнения принято делить на две группы: обыкновенные дифференциальные уравнения (ОДУ) и уравнения в частных производных. В данной лекции рассматриваются методы решения задач, описываемых обыкновенными дифференциальными уравнениями. Эти уравнения содержат только одну независимую переменную, в качестве которой может выступать время или пространственная координата. Иначе говоря, в таких уравнениях все функции зависят только от одной переменной и их производные по этой переменной являются полными. Уравнения в частных производных содержат более одной независимой переменной. Этими переменными могут быть, например, одновременно пространственные координаты и время или только пространственные координаты для статической задачи. В таких уравнениях производные от функций по любой из независимых переменных являются частными. Кроме того, уравнение может содержать смешанные производные. 10.3. Задача Коши Важным элементом задач, содержащих дифференциальные уравнения, являются дополнительные условия, которые необходимы для получения количественного решения. Применительно к обыкновенным дифференциальным уравнениям различают два вида задач: задачу с начальными условиями (задачу Коши) и задачу с краевыми условиями (краевую задачу). Задачу Коши можно сформулировать следующим образом. Дано обыкновенное дифференциальное уравнение y′ = dy(t)/dt = f (y, t) (10.1) и начальное условие y0 = y(t0) . (10.2) Требуется найти функцию y(t), удовлетворяющую уравнению (10.1) и начальному условию (10.2). На практике подобные задачи обычно связаны с расчѐтом переходных электрических, нестационарных тепловых или механических процессов при заданном в некоторый начальный момент времени исходном состоянии системы. 10.4. Сущность процедуры численного решения ОДУ Решение дифференциального уравнения (10.1) численным методом сводится к следующему: для заданной последовательности значений аргумента t0, t1, t2, …, tn и начального условия y0 = y(t0) определяется набор значений функции y0, y1, y2, …, yn (без определения ее аналитического вида y = y(t))), удовлетворяющих условиям: y0 = y(t0) , yk = y(tk) , k = 1, 2, …, n . (10.3) Рассмотрим три наиболее распространенных при решении практических задач численных метода интегрирования Эйлера, РунгеКутта и Адамса. Все они базируются на конечно-разностных методах представления производных. 81 10.5. Одношаговые конечно-разностные методы решения задачи Коши К одношаговым конечно-разностным методам относятся метод Эйлера (первого порядка), его модификация (второго порядка) и методы Рунге–Кутта (более высоких порядков). Понятие конечно-разностных методов представления производных было рассмотрено в предыдущих лекциях, посвященных методам интерполяции и численному вычислению производных. Напомним принципы, заложенные в этих методах применительно к дифференциальным уравнениям. Конечно-разностными называют уравнения, в которых бесконечно-малые приращения изменяющихся величин (dt) заменены достаточно малыми конечными приращениями (t). Большинство наиболее распространенных численных методов решения дифференциальных уравнений реализованы в виде рекуррентных формул, в которых очередной шаг решения осуществляется с использованием данных, полученных на одном или нескольких предыдущих шагах. Их принято подразделять на две большие группы:  Одношаговые методы, в которых для нахождения следующей точки на кривой, используется информация лишь об одном предыдущем шаге. К данной группе относятся методы Эйлера и Рунге-Кутта.  Многошаговые (методы прогноза и коррекции), в которых для отыскания следующей точки на кривой, требуется информация более чем об одной из предыдущих точек. В них для получения достаточно точного решения часто прибегают к итерации. К числу таких методов относятся методы Милна, Адамса-Бошфорта, Хемминга и ряд других. 10.5.1. Метод Эйлера Этот метод является простейшим численным методом решения задачи Коши. Рассмотрим его на примере решения обыкновенного дифференциального уравнения первого порядка (10.1) с соответствующим начальным условием (10.2). Заменим приращение в выражении (10.1) дифференциал dt на малое, но конечное приращение t = h. Тогда соответствующее ему приращение y будет равно: y =h · f (t, y). Используя начальное условие (10.2), получаем новое значение y: y1 = y0 +y = y0 + h · f (t0 , y0). Распространяя этот подход на последующие шаги решения, получим конечноразностную формулу численного решения задачи Коши в виде: yi+1 = yi + h · f (ti , yi). (10.4) Эта формула известна как формула метода Эйлера (часто ее называют явным методом Эйлера). Здесь (рис. 10.1) положение новой точки определяется по наклону кривой, вычисленному с помощью производной в данной точке. Таким образом, график численного решения представляет собой последовательность 82 коротких прямолинейных отрезков, которыми аппроксимируется истинная кривая y = y(x). Сам численный метод определяет порядок действий при переходе от данной точки к следующей. Рис. 10.1. Схема решения ОДУ явным методом Эйлера Отметим, что данную формулу можно получить, используя разложение функции y(t) в ряд Тейлора в окрестности некоторой точки ti с последующим отбрасыванием всех членов с дифференциалами высших порядков за исключением первого. Соответственно, использование данного подхода приводит к ошибке аппроксимации, порядок которой соответствует Rn(h2). Несмотря на достаточно большую погрешность данного метода, его простота в сочетании с физической и математической прозрачностью обусловили его широкое применение на практике. Кроме того на базе этого метода легче понять алгоритмы функционирования более сложных методов. В настоящее время этот метод используется при решении ряда физических задач, для которых погрешность порядка долей процента и даже нескольких процентов в ряде случаев является приемлемой. Это позволяет решать данным методом многие физические задачи достаточно просто и наглядно, попросту выбирая приемлемый для данной задачи шаг h. 10.5.2. Модифицированный метод Эйлера В рассмотренном явном методе Эйлера (10.4) для снижения погрешности расчетов требуется уменьшение величины шага интегрирования h, что приводит к увеличению процесса решения задачи (иначе времени интегрирования). Хотя тангенс наклона касательной к истинной кривой в исходной точке известен и равен y′(t0), он изменяется в соответствии с изменением независимой переменной. Поэтому в точке t0 + h наклон касательной уже те таков, каким он был в точке t0. Следовательно, при сохранении начального наклона касательной во всем интервале h в результаты вычислений вносится погрешность (рис. 10.1). Точность метода Эйлера можно существенно повысить, улучшив алгоритм аппроксимации 83 производной. Это можно сделать, например, используя среднее значение производной в начале и конце интервала. Именно это положено в основу модифицированного метода Эйлера (неявного метода Эйлера), алгоритм которого включает выполнение следующих шагов (рис. 10.2): 1. Сначала вычисляется значение функции в следующей точке по явному методу Эйлера: yn+1 = yn + h·f(tn, yn) , 2. Полученное значение используется для приближенного вычисления производной в конце интервала f(tn+1 , yn+1 ). 3. Вычисляется среднее значение производной между ее значениями в начале и конце интервала с последующей подстановкой ее в обычную формулу Эйлера (10.4) для получения более точного значения yn+1 yn+1 = yn + h·[f(tn, yn) + f(tn+1 , yn+1 )]/2 . (10.5) Рис. 10.2. Схема решения ОДУ модифицированным методом Эйлера Принцип, на котором основан модифицированный метод Эйлера, можно пояснить иначе. Для этого вернемся к разложению функции в ряд Тейлора Очевидно, что повышение точности расчета может быть достигнуто за счет сохранения члена, содержащего h2, и отбрасыванием членов более высоких порядков. Однако, чтобы сохранить член с h2 , надо знать вторую производную y(t0) . Ее можно аппроксимировать конечной разностью: Подставив это выражение в ряд Тейлора с отброшенными членами более высоких порядков, найдем (10.6) 84 что совпадает с выражением (10.5) (при условии, что точное значение производной y’(x + h) заменяется приближенным f(xn+1, y*n+1). Оценки показывают, что ошибка при этом имеет порядок h3). Этот метод является методом второго порядка, т.к. в нем используется член ряда Тейлора, содержащий h2. Ошибка на каждом шаге при использовании этого метода имеет порядок h3. За повышение точности приходится расплачиваться дополнительными затратами машинного времени, необходимого для вычисления yn+1 . Более высокая точность может быть достигнута при сохранении большего числа членов ряда Тейлора. Эта идея лежит в основе методов Рунге-Кутта. 10.5.3. Метод Рунге-Кутта Данный метод является одним из наиболее распространенных численных методов интегрирования обыкновенных дифференциальных уравнений. По сравнению с описанным выше методами Эйлера метод Рунге-Кутта имеет более высокую точность. Пусть на отрезке [a, b] требуется найти численное решение задачи Коши (10.1), где a = t0. Как и в предыдущем методе разобьем этот участок на n равных частей и построим последовательность значений t0, t1, t2, …, tn аргумента t искомой функции y(t). Предполагаем существование непрерывных производных функции y(t) до пятого порядка. Формулу для решения можно записать в виде: yk+1 = yk + yk (10.7) где yk — приращение искомой функции y(t) на (k+1)-ом шаге интегрирования. Придадим аргументу t приращение, равное шагу интегрирования h, и разложим функцию y(t + h) в ряд Тейлора в окрестности точки t, сохранив в нем пять членов: . Перенося первое слагаемое в этой сумме в левую часть получим, что (10.8) Здесь производные дифференцированием уравнения y = y (t) . определяются последовательным Вместо непосредственных вычислений по формуле (10.6) в методе Рунге-Кутта для каждого значения yk = y(xk ) определяются четыре числа: (10.9) 85 Доказывается, что если числа k1k, k2k, k3k, k4k, последовательно умножить на 1/6, 1/3, 1/3, 1/6 и сложить между собой, то выражение (10.10), соответствующее формуле Рунге-Кутта: 1 y k  (k1k  2k 2k  2k3k  k 4k ) 6 . (10.10) имеет погрешность Rn(h5). Таким образом, рабочая формула Рунге-Кутта для интегрирования имеет вид: (10.11) В отличие от расчетной схемы метода Эйлера, в которой каждое следующее значение yk+1 вычисляется непосредственно по единой формуле, в методе РунгеКутта необходимо проведение промежуточных вычислений по формулам (10.9) и (10.10). Примечание. Метод Рунге-Кутта может быть использован и при решении систем дифференциальных уравнений. Рассмотрим задачу Коши для системы второго порядка). В этом случае приращения y k и z k вычисляются по формулам: 1 y k  (k1k  2k 2 k  2k 3k  k 4k ) 6 1 z k  (m1k  2m2k  2m3k  m4k ) 6 (10.12) где k1k  hf1 ( xk , y k , z k ); m1k  hf 2 ( xk , y k , z k ); k 2 k  hf1 ( xk  0,5h , y k  0,5k1k , z k  0,5m1k ); m2 k  hf 2 ( xk  0,5h , y k  0,5k1k , z k  0,5m1k ); k 3k  hf1 ( xk  0,5h , y k  0,5k 2 k , z k  0,5m2 k ); (10.13) m3k  hf 2 ( xk  0,5h , y k  0,5k 2 k , z k  0,5m2 k ); k 4 k  hf1 ( xk  h , y k  k 3k , z k  m3k ); m4 k  hf 2 ( xk  h , y k  k 3k , z k  m3k ). Приближенное интегрирование системы уравнений (22) осуществляется по формулам вида: y k 1  y k  y k , z k 1  z k  z k , k  0, n (10.14) 10.5.4. Метод Адамса Пусть для задачи Коши найдены каким-либо способом (например, методом Эйлера или Рунге-Кутта) три последовательных значения искомой функции y1  y( x1 )  y( x0  h); y2  y( x2 )  y( x0  2h); y3  y( x3 )  y( x0  3h); Вычислим величины 86 q0  hy0'  hf ( x0 , y0 ) , q1  hy1'  hf ( x1 , y1 ) , q2  hy2'  hf ( x2 , y 2 ) , q3  hy3'  hf ( x3 , y3 ) . Метод Адамса позволяет найти решение задачи – функцию y(t) – в виде таблицы функций. Продолжение полученной таблицы из четырех точек осуществляется по экстраполяционной формуле Адамса: 1 5 3 y k 1  y k  qk  qk 1  2 qk 2  3 qk 3 , k  3,4,... 2 12 8 Затем уточнение проводится по интерполяционной формуле Адамса: 1 1 1 y k 1  y k  qk  qk 1  2 qk 2  3 qk 3 , k  3,4,.... (10.15) 2 12 24 Метод Адамса легко распространяется и на системы дифференциальных уравнений. Погрешность метода Адамса имеет тот же порядок, что и для метода Рунге-Кутта. КОНТРОЛЬНЫЕ ВОПРОСЫ Вопрос 1. Сформулируйте задачу Коши. Вопрос 2. Поясните сущность процедуры численного решения ОДУ. Вопрос 3. Перечислите рассмотренные в лекции одношаговые методы численного интегрирования. Вопрос 4. Сравните между собой явные и неявные методы численного интегрирования. Вопрос 5. Объясните процедуру методов численного решения ОДУ второго и более высокого порядка на примере метода Рунге-Кутта. 87 Лекции 11, 12 ОСОБЕННОСТИ И СПЕЦИФИКА ПРЕДМЕТА ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ Содержание  Исторический экскурс  Предмет вычислительной математики  Обусловленность задачи  Влияние выбора вычислительного алгоритма на результаты вычислений  Экономичность вычислительного метода  Погрешность метода  Элементы теории погрешностей  Задача численного дифференцирования 11.1. Исторический экскурс Первое применение вычислительных методов принадлежит древним египтянам, которые умели вычислять диагональ квадрата за конечное количество действий. Они также могли находить квадратный корень из 2, скорее всего, с помощью алгоритма, в дальнейшем получившего название формулы Герона, а еще позднее — метода Ньютона: С именем среднеазиатского врача, философа и математика Аль-Хоремзи связано понятие алгоритма. Разработкой вычислительных методов занимались Л.Эйлер, которому принадлежит, по-видимому, первый численный метод для решения обыкновенных дифференциальных уравнений, И.Ньютон, О.Л.Коши, Ж.Л.Лагранж, А.М.Лежандр, П.С.Лаплас, А.Пуанкаре, П.Л.Чебышѐв и многие другие известные математики. Решающую роль в развитии вычислительной математики как самостоятельной науки сыграли немецкий математик Карл Рунге и русский математик, механик и кораблестроитель А.Н.Крылов. В наше время вычислительная математика получила значительный импульс в 1950е годы, что было связано с развитием ядерной физики, механики полета, аэродинамики спускаемых космических аппаратов. В дальнейшем решались задачи, связанные не только с расчетами действия ядерного взрыва и обтеканием боеголовок стратегических ракет. Численные методы нашли свое применение в таких областях как динамика атмосферы, термогидрография, физика плазмы, механика горных пород и ледников, синергетика, биомеханика, теория оптимизации, математическая экономика и др. Наиболее наукоемки и требуют максимальных вычислительных ресурсов задачи физики, механики и электродинамики сплошных сред. К ним относятся системы уравнений в частных производных Эйлера, Лагранжа, Максвелла и др., кинетической теории газов 88 (уравнения Власова, Берда), а также задачи многомерной оптимизации. Развитие многих вычислительных методов — заслуга ученых, неразрывно связанных с МФТИ. Среди них следует упомянуть академиков А.А. Дородницына, О.М.Белоцерковского, А.А.Самарского, членов-корреспондентов РАН А.С.Холодова, Б.Н.Четверушкина, Ю.П.Попова, С.П.Курдюмова, профессоров В.С.Рябенького, Э.Э.Шноля, Р.П.Федоренко, Л.А.Чудова, В.Ф.Дьяченко, И.М.Гельфанда, Г.А.Тирского, А.П.Фаворского. 11.2. Предмет вычислительной математики Вычислительная математика отличается от других математических дисциплин и обладает специфическими особенностями. 1. Вычислительная математика имеет дело не только с непрерывными, но и с дискретными объектами. Так, вместо отрезка прямой часто рассматривается система точек , вместо непрерывной функции f(x) — табличная функция , вместо первой производной — ее разностная аппроксимация, например, Такие замены, естественно, порождают погрешности метода. 2. В машинных вычислениях присутствуют числа с ограниченным количеством знаков после запятой из-за конечности длины мантиссы при представлении действительного числа в памяти ЭВМ. Другими словами, в вычислениях присутствует машинная погрешность (округления) . Это приводит к вычислительным эффектам, неизвестным, например, в классической теории обыкновенных дифференциальных уравнений, уравнений математической физики или в математическом анализе. 3. В вычислительной практике большое значение имеет обусловленность задачи, т.е. чувствительность ее решения к малым изменениям входных данных. 4. В отличие от "классической" математики выбор вычислительного алгоритма влияет на результаты вычислений. 5. Существенная черта численного метода — экономичность вычислительного алгоритма, т.е. минимизация числа элементарных операций при выполнении его на ЭВМ. 6. Погрешности при численном решении задач делятся на две категории — неустранимые и устранимые. К первым относят погрешности, связанные с построением математической модели объекта и приближенным заданием входных данных, ко вторым — погрешности метода решения задачи и ошибки округления, которые являются источниками малых возмущений, вносимых в решение задачи. 89 Специфику машинных вычислений можно пояснить на нескольких элементарных примерах. 11.3. Обусловленность задачи Пример 11.1. Вычислить все корни уравнения Точное решение задачи легко найти: (x - 2)2 = ± 10- 4, x1= 2,01; x2= 1,99; x3,4= 2 ± 0,01i. Если компьютер работает при , то свободный член в исходном уравнении будет округлен до 16,0 и, с точки зрения представления чисел с плавающей точкой, будет решаться уравнение (x-2)4= 0, т.е. x1,2,3,4 = 2, что, очевидно, неверно. В данном случае малые погрешности в задании свободного члена ≈10-8 привели, независимо от метода решения, к погрешности в решении ≈ 10-2. Пример 11.2. Решается задача Коши для обыкновенного дифференциального уравнения 2-го порядка: u''(t) = u(t), u(0) = 1, u'(0) = - 1. Общее решение имеет вид u(t) = 0,5[u(0) + u'(0)]et + 0,5[u(0) - u'(0)]e- t. При заданных начальных данных точное решение задачи: u(x) = e-t, однако малая погрешность в их задании приведет к появлению члена , который при больших значениях аргумента может существенно исказить решение. Пример 11.3. Пусть необходимо найти решение обыкновенного дифференциального уравнения Его решение: приближенно: , однако значение u(t0) известно лишь , и на самом деле . Соответственно, разность u* - u будет 90 Предположим, что необходимо гарантировать некоторую заданную точность вычислений ε > 0 всюду на отрезке условие . Тогда должно выполняться Очевидно, что Отсюда можно получить требования к точности задания начальных данных при t0= 0. Таким образом, требование к заданию точности начальных данных оказываются в e10 раз выше необходимой точности результата решения задачи. Это требование, скорее всего, окажется нереальным. Решение оказывается очень чувствительным к заданию начальных данных. Такого рода задачи называются плохо обусловленными. Пример 11.4. Решением системы линейных алгебраических уравнений (СЛАУ) является пара чисел {1, 1}. Изменив правую часть системы на 0,01, получим возмущенную систему с решением {11.01; 0.00}, сильно отличающимся от решения невозмущенной системы. Эта система также плохо обусловлена. Пример 11.5. Рассмотрим полином (x - 1)(x - 2)...(x - 20)=x20 - 210x19 + ..., корни которого x1 = 1, x2 = 2, …, x20 = 20. Положим, что коэффициент (-210) при x19 увеличен на ≈ 10-7. В результате вычислений с 11-ю значащими цифрами получим совершенно иные корни: x1 = 1,00; x2 = 2,00; x3 = 3,00; x4 = 4,00; x5 = 5,00; x6 = 6,00; x7 = 7,00; x8 = 8,01; x9 = 8,92; x10,11 = 10,1 ± 0,644i; x12,13 91 = 11,8 ± 1,65i; x14,15 = 14,0 ± 2,52i; x16,17 = 16,7 ± 2,81i; x18,19 = 19,5 ± 19,4i; x20 = 20,8. Причина значительного расхождения также заключается в плохой обусловленности задачи вычисления корней рассматриваемого выражения. 11.4. Влияние выбора вычислительного алгоритма на результаты вычислений Пример 11.6. Пусть необходимо вычислить значение выражения Избавившись от знаменателя, получаем Полагая а) в) и рассматривая эти приближения как разные методы вычисления, получим следующие результаты: 7/5 0,004096 0,008000 1 17/12 0,005233 0,004630 - 0,1(6) Очевидно, что столь значительное различие в результатах вызвано влиянием ошибки округления в задании . Пример 11.7. Вычисление функции sin x с помощью ряда Тейлора. Из курса математического анализа известно, что функция синус представляется своим рядом Тейлора причем радиус сходимости ряда равен бесконечности — ряд сходится при любых значениях x. 92 Вычислим значения синуса при двух значениях аргумента. Пусть сначала . Будем учитывать лишь члены ряда, большие, чем 10- 4. Выполнив вычисления с четырьмя значащими цифрами, получим sin (0.5236) = 0.5000, что соответствует принятой точности. Пусть теперь . Если вычисления по данной формуле проводить с восемью значащими цифрами, то получим абсурдный результат: sin (25,66) ≈ 24 (учитывались члены ряда, большие, чем 10-8). Разумеется, выходом из создавшейся ситуации может быть использование формул приведения. Пример 11.8. Вычисление функции ex с помощью ряда Тейлора. Из курса математического анализа известно, что экспонента представляется своим рядом Тейлора радиус сходимости этого ряда также равен бесконечности. Приведем некоторые результаты расчетов ( вычисленные на компьютере). x ex 1 2,718282 20 -1 0,3678795 -10 -20 — значения экспоненты, exM 2,718282 0,3678794 1,202966 Выходом из этой ситуации может быть использование для отрицательных аргументов экспоненты формулы Естественно ожидать рост ошибок округления при вычислении рассматриваемой функции при больших значениях аргумента x. В этом случае можно использовать формулу ex = en + a = enea, где n = [x]. Пример 11.9. Рассмотрим следующий метод вычисления интеграла Интегрирование по частям дает 93 откуда следует . Тогда Очевидно, что отрицательные значения при n = 9,10 не имеют смысла. Дело в том, что ошибка, сделанная при округлении I0 до 6-ти значащих цифр сохранилась при вычислении I1, умножилась на 2! при вычислении I2, на 3! — при вычислении I3, и так далее, т.е. ошибка растет очень быстро, пропорционально n!. Пример 1.10. Рассмотрим методический пример вычислений на модельном компьютере, обеспечивающем точность Проанализируем причину происхождения ошибки, например, при вычитании двух чисел, взятых с точностью до третьей цифры после десятичной точки u = 1,001, v = 1,002, разность которых составляет Δ = |vM - uM| = 0,001. В памяти машины эти же числа представляются в виде Тогда Относительная ошибка при вычислении разности uM - vM будет равна Очевидно, что могут оказаться неверными. т.е. все значащие цифры Пример 1.11. Рассмотрим рекуррентное соотношение ui+1 = qui, q > 0, ui > 0. , u0 = a, 94 Пусть при выполнении реальных вычислений с конечной длиной мантиссы на i-м шаге возникла погрешность округления, и вычисления проводятся с возмущенным значением , тогда вместо ui+1 получим , т.е. Следовательно, если |q| > 1, то в процессе вычислений погрешность, связанная с возникшей ошибкой округления, будет возрастать (алгоритм неустойчив). В случае погрешность не возрастает и численный алгоритм устойчив. 11.5. Экономичность вычислительного метода Пример 1.12. Пусть требуется вычислить сумму S = 1 + x + x2 + … + x1023 при 0 < x < 1. Для последовательного вычисления необходимо проделать 1022 умножения, а затем столько же сложений. Однако если заметить, что то количество арифметических действий значительно уменьшается; в частности, для вычисления x1024 требуется всего 10 умножений: Пример 1.13. Вычисления значений многочленов. Если вычислять значение многочлена P(x) = a0 + a1x + a2x2 + … + anxn "в лоб", т.е. вычислять значения каждого члена и суммировать, то окажется, что необходимо выполнить (n2 + [n/2]) умножений и n сложений. Кроме того, такой способ вычислений может привести к накоплению ошибок округления при вычислениях с плавающей точкой. Его очевидным улучшением является вычисление каждого члена последовательным умножением на x. Такой алгоритм требует (2n - 1) умножение и n сложений. Еще более экономичным алгоритмом является хорошо известная в алгебре схема Горнера: P(x) = ((...((anx + an - 1)x + an - 2)x + ... + a0), требующая n операций сложения и n операций умножения. Этот метод был известен в средние века в Китае под названием Тянь-Юань и был заново открыт в Европе в начале XIX века англичанином Горнером и итальянцем Руффини. 95 Пример 1.14. Рассмотрим систему линейных алгебраических уравнений (СЛАУ) вида матрицей , с трехдиагональной Если проводить решение такой системы без учета специфической структуры матрицы (например, с помощью метода Гаусса), то количество арифметических действий будет порядка n3, если же учесть эту структуру, то количество операций можно уменьшить до n. 11.6. Погрешность метода Оценим погрешность при вычислении первой производной при помощи соотношения: где O(h) есть погрешность метода. В данном случае под погрешностью метода понимается абсолютная величина разности которая составляет O(h) (более точно где ). Если же взять другой метод вычисления производной то получим, что его погрешность составляет O(h2), это оказывается существенным при малых h. Однако уменьшать h до бесконечности не имеет смысла, что видно из следующего примера. Реальная погрешность при вычислении первой производной будет поскольку абсолютная погрешность вычисления значения функции за счет машинного округления не превосходит 96 В этом случае можно найти оптимальный шаг h. Будем считать полную погрешность в вычислении производной Δ функцией шага h. Отыщем минимум этой функции. Приравняв производную Δ'h(h) к нулю, получим оптимальный шаг численного дифференцирования Выбирать значение h меньше оптимального не имеет смысла, так как при дальнейшем уменьшении шага суммарная погрешность начинает расти из-за возрастания вклада ошибок округления. 11.7. Элементы теории погрешностей Определение. Пусть u и u* — точное и приближенное значение некоторой величины, соответственно. Тогда абсолютной погрешностью приближения u* называется величина Δ (u*), удовлетворяющая неравенству Определение. Относительной погрешностью называется величина удовлетворяющая неравенству , Обычно используется запись Определение. Пусть искомая величина u является функцией параметров , u* — приближенное значение u. Тогда предельной абсолютной погрешностью называется величина Предельной относительной погрешностью называется величина D(u*)/| u*|. Пусть — приближенное значение Предполагаем, что u - непрерывно дифференцируемая функция своих аргументов. Тогда, по формуле Лагранжа, где 97 Отсюда где . Можно показать, что при малых эта оценка не может быть существенно улучшена. На практике иногда пользуются грубой (линейной) оценкой Несложно показать, что a) , предельная погрешность суммы или разности равна сумме предельных погрешностей. b) Предельная относительная погрешность произведения или частного приближенного равна сумме предельных относительных погрешностей 11.8. Задача численного дифференцирования В пункте 1.2 уже была введена простейшая формула численного дифференцирования. Рассмотрим задачу приближенного вычисления значения производной подробнее. Пусть задана таблица значений xi. В дальнейшем совокупность точек на отрезке, на котором проводятся вычисления, иногда будут называться сеткой, каждое значение xi — узлом сетки. Пусть сетка — равномерная, и расстояние между узлами равно h — шагу сетки. Пусть узлы сетки пронумерованы в порядке возрастания, т.е. x0 = a, xj = a + jh,, j = 0, 1, ... , N Пусть f(xj) = fj — функция, определенная в узлах сетки. Такие функции будут называться табличными, или сеточными функциями. Считаем, кроме того, что рассматриваемая сеточная функция есть проекция (или ограничение) на сетку некоторой гладкой нужное число раз непрерывно дифференцируемой функции f(x). По определению производной 98 тогда, если шаг сетки достаточно мал, по аналогии можно написать формулу в конечных разностях, дающую приближенное значение производной сеточной функции: (1.1) Если параметр h достаточно мал, то можно считать полученное значение производной достаточно точным. Погрешность формулы (1.1) оценена в пункте 1.2. Как показано выше, при уменьшении шага сетки h ошибка будет уменьшаться, но при некотором значении h ошибка может возрасти до бесконечности. При оценке погрешности метода обычно считается, что все вычисления были точными. Но существует ошибка округления. При оценке ее большую роль играет машинный ε — мера относительной погрешности машинного округления, возникающей из-за конечной разрядности мантиссы при работе с числами в формате с плавающей точкой. Напомним, что по определению машинным ε называют наибольшее из чисел, для которых в рамках используемой системы вычислений выполнено 1 + ε = 1. Тогда абсолютная погрешность при вычислении значения функции (или представлении табличной функции) есть . Максимальный вклад погрешностей округления при вычислении производной по формуле (1.1) будет тогда, когда члены в знаменателе (1.1) имеют ошибки разных знаков. Пусть k = max |f(x)|, максимум ищется на отрезке, на котором вычисляются значения производных. Тогда суммарная ошибка, состоящая из погрешности метода и погрешности округления, есть Для вычисления оптимального шага численного дифференцирования найдем минимум суммарной ошибки, как функции шага сетки откуда Если требуется повысить точность вычисления производных, необходимо воспользоваться формулами, имеющими меньшие погрешности метода. Так, из курсов математического анализа известно, что По аналогии напишем конечно-разностную формулу 99 (1.2) (1.2) — формула с центральной разностью. Исследуем ее на аппроксимацию, т.е. оценим погрешность метода. Предположим, что функция, которую спроектировали на сетку, трижды непрерывно дифференцируема, тогда Погрешность метода определяется 3-й производной функции. Введем , тогда суммарная погрешность при вычислении по формуле с центральной разностью есть для вычисления оптимального шага, находя минимум погрешности, как функции шага сетки, имеем откуда Для более точного вычисления производной необходимо использовать разложение более высокого порядка, шаг hopt будет увеличиваться. Формула (1.1) — двухточечная, (1.2) — трехточечная: при вычислении производной используются точки (узлы) xj (узел входит с нулевым коэффициентом), xj + h, xj - h — совокупность узлов, участвующих в каждом вычислении производной, в дальнейшем будем иногда называть сеточным шаблоном. Введем на рассматриваемом отрезке шаблон из нескольких точек. Считаем, что сетка равномерная — шаг сетки постоянный, расстояния между любыми двумя соседними узлами равны. Используем для вычисления значения первой производной следующую приближенную (конечно-разностную) формулу: (1.3) шаблон включает l точек слева от рассматриваемой точки xj и m справа. Коэффициенты — неопределенные коэффициенты. Формула дифференцирования может быть и односторонней — либо l, либо m могут равняться нулю. В первом случае иногда называют (на наш взгляд, не слишком 100 удачно) такую приближенную формулу формулой дифференцирования вперед, во втором — формулой дифференцирования назад. Потребуем, чтобы (1.3) приближала первую производную с точностью O(hl + m) . Используем разложения в ряд Тейлора в окрестности точки xj. Подставляя их в (1.3), получим Потребуем выполнение условий: (1.4) Получаем систему линейных алгебраических уравнений для неопределенных коэффициентов (1.4). Матрица этой системы есть Вектор правых частей (0, 1, 0, ..., 0)T. Определитель данной матрицы — детерминант Вандермонда. Из курса линейной алгебры следует, что он не равен нулю. Тогда существует единственный набор коэффициентов , который позволяет найти на шаблоне из (1 + l + m) точек значение первой производной с точностью O(hl + m). Для нахождения второй производной можно использовать ту же самую формулу (1.3) с небольшой модификацией только теперь Очевидно, что и данная система уравнений для нахождения неопределенных коэффициентов имеет единственное решение. Для получения с той же точностью приближенных значений производных до порядка l + m включительно с точностью O(hl + m) модификации формулы (1.3) и условий (1.4) очевидны, набор неопределенных коэффициентов находится единственным образом. 101 Таким образом, доказано следующее утверждение. На сеточном шаблоне, включающем в себя N + 1 точку, с помощью метода неопределенных коэффициентов всегда можно построить единственную формулу для вычисления производной от первого до n порядка включительно с точностью O(hN). Утверждение доказано для равномерной сетки, но на случай произвольных расстояний между сеточными узлами обобщение проводится легко. Так как на практике вычисления проводятся с конечной длиной мантиссы, то получить нулевую ошибку невозможно. Читателям предлагается оценить значение оптимального шага при вычислениях по формулам типа (1.3) самостоятельно. 1.7. Задачи 1. Найти абсолютную предельную погрешность, погрешность по производной, линейную погрешность для функции u = t10, если заданы точка приближения t* = 1, значение функции u* в этой точке и погрешность Δt* = 10 -1. Решение : обозначим Абсолютная предельная погрешность может быть определена как Оценка погрешности u при вычислении значения функции по максимуму производной и линейная оценка соответственно будут D1(u*) = bΔ(t*) = 2,3 ...; D2(u*) = |γ(0)|Δ(t*) = 1. 2. Дать линейную оценку погрешности при вычислении неявной функции , если известны точка приближения , * значение функции в точке приближения u и погрешность в определении аргументов . Решение. Дифференцируя по tj, получим откуда 102 При заданных , можно найти u* как корень уравнения , а затем — значения откуда можно получить линейную оценку погрешности функции D2(u*). 3. Вычислить относительную погрешность в определении значения функции 4. u = xy2z3, если x* = 37,1, y* = 9,87, z* = 6,052, Δx* = 0,3, Δy* = 0,11, 5. Δz* = 0,016. Решение: 6. Оценить погрешность в определении корней квадратного уравнения , если заданы приближения . Пусть u* — решение уравнения Из формулы получим Следовательно, линейная оценка будет 103 1.8. Задачи для самостоятельного решения 1. Найти абсолютную предельную погрешность, погрешность по производной и линейную оценку погрешности для функций u = sin t, Заданы точка приближения t = t* и погрешность Δ t. 2. Определить шаг , при котором погрешность вычисления производной u'(t), приближенно вычисляемой в соответствии с формулами не превосходит 103. Известно, что для любых t. 3. Пусть для вычисления функции u = f(t) используется частичная сумма ряда Маклорена причем аргумент задан с погрешностью Δ t = 10 -3 . Найти n такое, чтобы погрешность в определении функции u(y) по данной формуле не превышала Δ t. Рассмотреть отрезки . Предложить более совершенный алгоритм для вычисления функций u(t) = sin t, u(t) = et на отрезке . 4. Определить оптимальный шаг численного дифференцирования τ при использовании для вычисления производной приближенной формулы имеющей четвертый порядок точности, если известно, что значения функций вычисляются с точностью ε. ,а 5. Вычислить относительную погрешность в определении значения функции u(x,y,z) = x2y2/z4, если заданы 6. x* = 37,1, y* = 9,87, z* = 6,052, Δ(x*) = 0,1; 7. Δ(y*) = 0,05; Δ(z*) = 0,02. 104 КОНТРОЛЬНЫЕ ВОПРОСЫ Вопрос 1. Дайте краткое пояснение назначения и сфер применения предмета вычислительной математики. Вопрос 2. Объясните сущность понятия обусловленность задачи. Вопрос 3. Прокомментируйте связь выбора вычислительного алгоритма с результатом вычислений. Вопрос 4. Что понимается под экономичностью вычислительного метода? Вопрос 5. В чем заключается сложность задачи численного дифференцирования? 105