Главная
База знаний "Allbest"
Программирование, компьютеры и кибернетика
Методы компьютерных вычислений и их приложение к физическим задачам
Методы компьютерных вычислений и их приложение к физическим задачам
Точность вычислений, классификация погрешностей. Оценка апостериорной погрешности, численное дифференцирование. Численное решение систем линейных уравнений. Аппроксимация функций методом наименьших квадратов. Решение нелинейных уравнений с неизвестным.
| Рубрика | Программирование, компьютеры и кибернетика |
| Вид | методичка |
| Язык | русский |
| Дата добавления | 10.10.2010 |
| Размер файла | 611,8 K |
Отправить свою хорошую работу в базу знаний просто. Используйте форму, расположенную ниже
Студенты, аспиранты, молодые ученые, использующие базу знаний в своей учебе и работе, будут вам очень благодарны.
В. Пикулев
Методы компьютерных вычислений и их приложение к физическим задачам
Методическое пособие
Петрозаводск, 2004
1. Всякий специалист стремится достичь
своего уровня некомпетентности.
2. Компьютер - средство многократного
увеличения некомпетентности человека.
3. Если две ошибки не принесли результата
- испробуй третью.
Лоуренс Питер, «Иерархиология».
Введение в предмет
Численные методы - раздел математики, который со времен Ньютона и Эйлера до настоящего времени находит очень широкое применение в прикладной науке. Традиционно физика является основным источником задач построения математических моделей, описывающих явления окружающего мира, она же является основным потребителем алгоритмов и программ, позволяющих эти задачи с определенным успехом решать.
При этом задачей физика является не только правильный выбор программы, которая призвана решать физическую проблему, но и подробный анализ и корректировка используемых алгоритмов, в соответствии с реалиями поставленной задачи и теми математическими правилами, которые либо допускают существование решения с заданной точностью, либо говорят о невозможности такого решения.
Примеры современных физических задач, для решения которых используются численные методы - моделирование астрономических событий (рождение и развитие Вселенной), моделирование процессов в микромире (распад и синтез частиц), моделирование установок и процессов термоядерного синтеза. Более «прикладные» задачи - моделирование физических процессов в твердотельных структурах (широко используется в проектировании и изготовлении интегральных схем), моделирование процессов в газах и плазме. Учитывая большую сложность и дороговизну современных экспериментальных методик, и, с другой стороны, постоянный рост производительности вычислительных систем, нетрудно определить тенденцию к увеличению в настоящее время доли модельных (вычислительных) экспериментов. Большое количество численных методов разработано для решения задач математической физики, к которым, например, относятся задачи тепло- и массопереноса, исследования турбулентного движения.
К инженерным приложениям численных методов можно отнести расчеты магнитных и электростатических линз для заряженных частиц, различного рода радиотехнические расчеты, включая, например, проектирование СВЧ-волноводов. Любопытно, что как в теоретической физике, так и в инженерной практике решаются численными методами различные задачи теоретической механики, например, задачи столкновения (в том числе динамический хаос). Естественно, такое приложение вычислительной техники к физике, как управление экспериментом и сбор данных, в данном курсе не рассматривается.
План построения вычислительного эксперимента:
1) Создание модели, фиксирующей главные исследуемые факторы. Одновременно формулируются рамки применимости модели.
2) Предварительное исследование математической модели: поверка корректности постановки задачи, существования и единственности решения.
3) Разработка метода расчета сформулированной задачи, построение эффективных вычислительных алгоритмов.
4) Создание программы, осуществляющей моделирование физического объекта, включающей в себя реализации используемых численных методов, проверки корректности ввода исходных и вывода результирующих данных
5) Сравнение полученных результатов моделирования с тестовыми примерами и экспериментальными данными; решение вопроса о правильности практического моделирования (иначе повторяются пункты 3 и 4).
6) Решение вопроса о достоверности предложенной математической модели. Если модель не описывает экспериментальные данные, возврат на пункт 1.
Таким образом, численный эксперимент - это не однократное вычисление по некоторому набору формул, а многостадийный процесс программирования, анализа результатов и их погрешностей.
Задача называется корректно поставленной, если для любых входных данных х из некоторого класса решение y существует, единственно и устойчиво по входным данным.
(Класс может представлять собой координатное бесконечномерное пространство, множество непрерывных функций и др.)
Численный алгоритм -- однозначная последовательность действий, которые могут привести к одному решению.
Отсутствие устойчивости обычно означает, что сравнительно небольшой погрешности дx соответствует весьма большое дy, а значит получаемое решение будет далеко от истинного. К такой задаче численные методы применять бессмысленно, ибо погрешности численного расчета будут катастрофически нарастать. Устойчивость задачи определяется (1) математической формулировкой, (2) используемым алгоритмом расчета. Пример неустойчивой задачи в первом случае:
Система имеет решение , однако
Система имеет решение , то есть разница в коэффициенте менее 1% приводит к изменению решения в 300%.
Пример алгоритмической неустойчивости - вычисление производных численными методами: какой бы метод мы не использовали, приходится вычитать весьма мало различающиеся числа.
В настоящее время развиты методы решения многих некорректных задач, которые основаны на решении вспомогательной корректной задачи, близкой к исходной.
Если выполняется условие для норм (модулей), то задача устойчива. Однако если константа С очень велика, то фактически наблюдается слабая устойчивость. Такую задачу называют плохо обусловленной. Пример - дифференциальное уравнение с начальными условиями . Общее решение дифференциального уравнения есть . Начальные условия приводят к обнулению первого слагаемого, но если из-за погрешности начальных данных это будет не так, то при возрастании x влияние первого слагаемого будет катастрофически нарастать.
1 Точность вычислений, классификация погрешностей
Во всех случаях математическая точность решения должна быть в 2-4 раза выше, чем ожидаемая физическая точность модели. Более высокая математическая точность, как и более низкая, будут неадекватны данной модели.
Существуют четыре источника погрешности результата:
1) погрешность математической модели - связана с ее несоответствием физической реальности, так как абсолютная истина недостижима. Если математическая модель выбрана недостаточно тщательно, то, какие бы методы мы не применяли для расчета, все результаты будут недостаточно надежны, а в некоторых случаях и совершенно неправильны.
2) погрешность исходных данных, принятых для расчета. Это неустранимая погрешность, но это погрешность возможно и необходимо оценить для выбора алгоритма расчета и точности вычислений. Как известно, ошибки эксперимента условно делят на систематические, случайные и грубые, а идентификация таких ошибок возможно при статистическом анализа результатов эксперимента.
3) погрешность метода - основана на дискретном характере любого численного алгоритма. Это значит, что вместо точного решения исходной задачи метод находит решение другой задачи, близкого в каком-то смысле (например по норме банахова пространства) к искомому. Погрешность метода - основная характеристика любого численного алгоритма. Погрешность метода должна быть в 2-5 раз меньше неустранимой погрешности.
4) погрешность округления -- связана с использованием в вычислительных машинах чисел с конечной точностью представления.
Вот иллюстрация этих определений. Пусть имеется реальный маятник, совершающий затухающие колебания, начинающий движение в момент t = t0. Требуется найти угол отклонения ц от вертикали в момент t1. Движение маятника мы можем описать следующим дифференциальным уравнением:
,
где l - длина маятника, g - ускорение силы тяжести, м - коэффициент трения.
Как только принимается такое описание задачи, решение уже приобретает неустранимую погрешность, в частности потому, что реальное трение зависит от скорости не совсем линейно (погрешность модели). Кроме того, воспроизведя реальный эксперимент, мы зададим l, g (в известной точке планеты), м с некоторой точностью, и получим набор значений с погрешностью, которую можем оценить из анализа статистики некоторого числа однотипных опытов (погрешность исходных данных). Взятое в модели дифференциальное уравнение нельзя решить в явном виде, для его решения требуется применить какой-либо численный метод, имеющий заранее известную погрешность, которая должна быть меньше неустранимой погрешности. После совершения вычислений мы получим значения с погрешностью большей, нежели погрешность метода, так как к ней прибавится погрешность округления.
Рассмотрим правила расчета погрешности округления:
1) Сложение и вычитание приближенных чисел
Введем в рассмотрение два числа a и b, называемых приближенными, то есть это есть оценка точных значений A и B, известных с абсолютными погрешностями ±еa и ±еb. Знаки этих погрешностей нам неизвестны, следовательно для обеспечения достоверности конечного результата мы должны взять наихудший случай, когда погрешности складываются. Таким образом формулируются следующие правила:
1. Абсолютная погрешность суммы приближенных чисел равна сумме абсолютных погрешностей слагаемых.
2. Абсолютная погрешность разности приближенных чисел равна сумме абсолютных погрешностей слагаемых.
Относительной погрешностью приближенного числа a будет являться величина . По этому же правилу определим относительную погрешность суммы приближенных чисел a и b как . При этом можно показать, что
3. Относительная погрешность суммы слагаемых одного знака заключена между наименьшей и наибольшей относительными погрешностями слагаемых: .
4. Для разности двух приближенных чисел одного знака величина относительной погрешности может быть сколь угодно большой.
2) Умножение и деление приближенных чисел
Очевидно, что приближенное число . Тогда для произведения . Если пренебречь последним малым слагаемым в скобках, то можно сформулировать следующее правило:
1. Относительная погрешность произведения приближенных чисел равна сумме относительных погрешностей множителей .
Так как деление на число b равнозначно умножению на 1/b, то справедливо утверждение:
2. Относительная погрешность частного приближенных чисел равна сумме относительных погрешностей делимого и делителя.
Следовательно, при умножении и делении приближенных чисел необходимо принимать во внимание количество значащих цифр, характеризующих относительную точность числа, а не количество десятичных знаков, обуславливающих его абсолютную погрешность.
Совершенно очевидно, что при большом количестве действий такого сорта правила нельзя считать удовлетворительными, так как погрешности будут иметь разные знаки и компенсировать друг друга. Статистическая оценка показывает, что при N одинаковых действиях среднее значение суммарной ошибки больше единичной в раз, если нет систематических причин для накопления погрешности. Систематические причины возникают, если, например в алгоритме вычитаются близкие по величине числа.
При любых расчетах надо устанавливать такую точность вычислений, чтобы погрешность округления была существенно меньше всех остальных погрешностей.
2 Краткое введение в используемые программные средства
Традиционные языки высокоуровневого программирования. В большинстве практических случаев моделирование или численный расчет предполагает использование готового алгоритма, который необходимым образом модифицируется для конкретной задачи. В настоящее время существуют обширные фонды алгоритмов и программ, ориентированных на классические языки программирования, такие как Фортран, Си и Паскаль. Наиболее мощные математические библиотеки были разработаны для самого старого языка программирования из вышеперечисленных --- Фортрана. Однако во многих случаях, обладая соответствующей квалификацией, легче написать программу, отталкиваясь от первых принципов, нежели тестировать и исправлять, или переводить на другой язык программирования незнакомый текст.
Другое направление, которое в настоящее время особенно популярно - использование мощных математических пакетов для численных и аналитических расчетов.
Mathcad.
Разрабатывается компанией MathSoft Inc. Является наиболее легкой для освоения системой математических расчетов. Принята концепция «активного документа», то есть все вычисления записываются в традиционной математической нотации (с использованием значков интеграла, суммы и др.), а после введения знака равенства или другого запускающего символа появляется рассчитанное значение. Основной недостаток -- слишком мал набор основных функций и очень низкое быстродействие.
MathLab.
Система MATLAB (MATrix LABoratory) разрабатывается фирмой MathWorks. Эта система создана для работы в среде Windows и представляет собой интерактивную среду для вычислений и моделирования, причем она может работать как в режиме непосредственных вычислений (очень напоминает режим «командной строки»), так и в режиме интерпретации написанных программ. Сильная сторона системы - виртуозная работа с матрицами и векторами. Численное значение или аналитическая формула, а также сообщения системы выводится на экран в виде списка. Помимо обычных алгебраических вычислений система имеет огромный набор встроенных функций, а также имеется возможность создавать пользовательские функции. В системе очень качественно реализовано построение двух и трехмерных изображений, в том числе динамически изменяющихся. Кроме того, имеется библиотека, которая обеспечивает удобное управление исполнением программ. И это только базовый набор, который обычно расширяется многочисленными дополнениями - например языком Simulink моделирования нелинейных динамических систем. Основное назначение - технические расчеты.
Mathematica.
Разрабатывается фирмой Wolfram Research. Одна из наиболее сложных для освоения систем. Предназначена в основном для решения теоретических задач, в связи с чем весьма популярна в научных кругах. Предуставляет широкие возможности в проведении символических (аналитических) преобразований, красотой интерфейса не отличается.
Maple.
Разработка компании Waterloo Maple Inc. Также очень популярный в научных кругах пакет аналитических преобразований и численных методов. Символический процессор Maple поставляется отдельно, в связи с чем производители других программных средств могут интегрировать его в свои разработки.
3 Численное интегрирование
Задача численного интегрирования состоит в замене исходной подинтегральной функции f(x), для которой трудно или невозможно записать первообразную в аналитике, некоторой аппроксимирующей функцией ц(x). Такой функцией обычно является полином (кусочный полином) . То есть:
,
где - априорная погрешность метода на интервале интегрирования,
а r(x) - априорная погрешность метода на отдельном шаге интегрирования.
Обзор методов интегрирования.
Методы вычисления однократных интегралов называются квадратурными (для кратных интегралов - кубатурными).
1) Методы Ньютона-Котеса. Здесь ц(x) - полином различных степеней. Сюда относятся метод прямоугольников, трапеций, Симпсона.
2) Методы статистических испытаний (методы Монте-Карло). Здесь узлы сетки для квадратурного или кубатурного интегрирования выбираются с помощью датчика случайных чисел, ответ носит вероятностный характер. В основном применяются для вычисления кратных интегралов.
3) Сплайновые методы. Здесь ц(x) - кусочный полином с условиями связи между отдельными полиномами посредством системы коэффициентов.
4) Методы наивысшей алгебраической точности. Обеспечивают оптимальную расстановку узлов сетки интегрирования и выбор весовых коэффициентов с(x) в задаче. Сюда относится метод Гаусса-Кристоффеля (вычисление несобственных интегралов) и метод Маркова.
Метод прямоугольников.
Различают метод левых, правых и средних прямоугольников. Суть метода ясна из рисунка. На каждом шаге интегрирования функция аппроксимируется полиномом нулевой степени - отрезком, параллельным оси абсцисс.
Выведем формулу метода прямоугольников из анализа разложения функции f(x) в ряд Тейлора вблизи некоторой точки x = xi.
…
Рассмотрим диапазон интегрирования от xi до xi + h, где h - шаг интегрирования.
Вычислим …=
==. Получили формулу правых (или левых) прямоугольников и априорную оценку погрешности r на отдельном шаге интегрирования. Основной критерий, по которому судят о точности алгоритма - степень при величине шага в формуле априорной оценки погрешности.
В случае равного шага h на всем диапазоне интегрирования общая формула имеет вид
.
Здесь n - число разбиений интервала интегрирования, . Для справедливости существования этой оценки необходимо существование непрерывной f'(x).
Метод средних прямоугольников. Здесь на каждом интервале значение функции считается в точке , то есть . Разложение функции в ряд Тейлора показывает, что в случае средних прямоугольников точность метода существенно выше:
.
Метод трапеций.
Аппроксимация в этом методе осуществляется полиномом первой степени. Суть метода ясна из рисунка.
На единичном интервале.
В случае равномерной сетки (h = const)
При этом , а . Погрешность метода трапеций в два раза выше, чем у метода средних прямоугольников! Однако на практике найти среднее значение на элементарном интервале можно только у функций, заданных аналитически (а не таблично), поэтому использовать метод средних прямоугольников удается далеко не всегда. В силу разных знаков погрешности в формулах трапеций и средних прямоугольников истинное значение интеграла обычно лежит между двумя этими оценками.
Особенности поведения погрешности.
Казалось бы, зачем анализировать разные методы интегрирования, если мы можем достичь высокой точности, просто уменьшая величину шага интегрирования. Однако рассмотрим график поведения апостериорной погрешности R результатов численного расчета в зависимости от числа n разбиений интервала (то есть при шаг ). На участке (1) погрешность уменьшается в связи с уменьшением шага h. Но на участке (2) начинает доминировать вычислительная погрешность, накапливающаяся в результате многочисленных арифметических действий. Таким образом, для каждого метода существует своя Rmin, которая зависит от многих факторов, но прежде всего от априорного значения погрешности метода R.
Уточняющая формула Ромберга.
Метод Ромберга заключается в последовательном уточнении значения интеграла при кратном увеличении числа разбиений. В качестве базовой может быть взята формула трапеций с равномерным шагом h.
Обозначим интеграл с числом разбиений n = 1 как .
Уменьшив шаг в два раза, получим .
Если последовательно уменьшать шаг в 2n раз, получим рекуррентное соотношение для расчета .
Пусть мы вычислили четыре раза интеграл с n от 1 до 4. Представим следующий треугольник:
R(1;1)
R(2;1) R(2;2)
R(3;1) R(3;2) R(3;3)
R(4;1) R(4;2) R(4;3) R(4;4)
В первом столбце стоят значения интеграла, полученные при последовательном удвоении числа интервалов. Следующие столбцы - результаты уточнения значения интеграла по следующей рекуррентной формуле:
.
Правое нижнее значение в треугольнике - искомое уточненное значение интеграла.
Метод Симпсона.
Подинтегральная функция f(x) заменяется интерполяционным полиномом второй степени P(x) - параболой, проходящей через три узла, например, как показано на рисунке ((1) - функция, (2) -- полином).
Рассмотрим два шага интегрирования (h = const = xi+1 - xi), то есть три узла x0, x1, x2, через которые проведем параболу, воспользовавшись уравнением Ньютона:
.
Пусть z = x - x0,
тогда
Теперь, воспользовавшись полученным соотношением, сосчитаем интеграл по данному интервалу:
В итоге .
Для равномерной сетки и четного числа шагов n формула Симпсона принимает вид:
Здесь , а в предположении непрерывности четвертой производной подинтегральной функции.
4 Блок-схема алгоритма метода Симпсона
Методы Монте-Карло.
1)
1. одномерная случайная величина - статистический вариант метода прямоугольников.
В качестве текущего узла xi берется случайное число, равномерно распределенное на интервале интегрирования [a, b]. Проведя N вычислений, значение интеграла определим по следующей формуле:
.
Для R можно утверждать хотя бы ~.
2) двумерная случайная величина -- оценка площадей.
Рассматриваются две равномерно распределенных случайных величины xi и yi, которые можно рассматривать как координаты точки в двумерном пространстве. За приближенное значение интеграла принимается количества точек S, попавших под кривую y = f(x), к общему числу испытаний N, т.е. .
И первый, и второй случай легко обобщаются на кратные интегралы.
5 Оценка апостериорной погрешности
Мы записывали априорные оценки главного члена погрешности в виде
R0 = Ahp, (1)
где A - коэффициент, зависящий от метода интегрирования и вида подинтегральной функции; h - шаг интегрирования, p - порядок метода. Такого сорта оценку можно применить не только к методам интегрирования, но и ко многим другим численным алгоритмам.
Первая формула Рунге.
Пусть w - точное значение, к которому должен прийти численный метод (мы его не знаем). Результат численного расчета дает нам величину wh такую, что
. (2)
Теперь вычислим ту же величину w с шагом kh, где константа k может быть как больше, так и меньше единицы. Коэффициент A будет одинаковый, так как вычисление осуществляется одним и тем же методом. Получаем
. (3)
Приравняем правые части выражений (2) и (3) и пренебрежем бесконечно малыми величинами одинакового порядка малости.
.
Отсюда, учитывая (1), получим
. (4)
Эта формула, выражающая апостериорную оценку главного члена погрешности величины w путем двойного просчета с разным шагом, носит название первой формулы Рунге. При уменьшении шага главный член погрешности будет стремиться к полной погрешности R.
Вторая формула Рунге.
Так как модуль и знак апостериорной погрешности из формулы (4) известны, можно уточнить искомое значение
Это вторая формула Рунге. Однако теперь погрешность wcorr не определена, известно лишь, что она по модулю меньше R0.
Алгоритм Эйткена.
Способ оценки погрешности для случая, когда порядок метода p неизвестен. Более того, алгоритм позволяет опытным путем определить и порядок метода. Для этого в третий раз вычислим значение величины w с шагом k2h:
. (5)
Приравняем правые части выражений (5) и (3): . Отсюда:
. Подставим сюда значение R0 из (4):
.
Из этой формулы определяем знаменатель для (4). Кроме того, определяем порядок . Для правильно реализованных алгоритмов методов априорных и апостериорных порядки должны получиться совпадающими. Программная реализация формул Рунге позволяет вычислить определенные интегралы с заданной точностью, когда выбор необходимого числа разбиений интервала интегрирования осуществляется автоматически. Пример - уже рассмотренная ранее формула Ромберга.
6 Численное дифференцирование
Методы численного дифференцирования применяются, если исходную функцию f(x) трудно или невозможно продифференцировать аналитически. Например, эта функция может быть задана таблично. Задача численного дифференцирования - выбрать легко вычисляемую функцию (обычно полином) , для которой приближенно полагают .
Численное дифференцирование - некорректная задача, так как отсутствует устойчивость решения. При численном дифференцировании приходится вычитать друг из друга близкие значения функции. Это приводит к уничтожению первых значащих цифр, т.е. к потере части достоверных знаков числа. А так как значения функции обычно известны с определенной погрешностью, то все значащие цифры могут быть потеряны. На графике кривая (1) соответствует уменьшению погрешности дифференцирования при уменьшении шага; кривая (2) представляет собой неограниченно возрастающий (осциллирующий) вклад неустранимой погрешности исходных данных - значений функции y(x). Критерий выхода за оптимальный шаг при его уменьшении -- «разболтка» решения: зависимость результатов вычислений становится нерегулярно зависящей от величины шага.
Пусть введена как интерполяционный многочлен Ньютона. В этом случае для произвольной неравномерной сетки:
, для i = 0,1…n-1, интерполяция полиномом первой степени.
, интерполяция полиномом второй степени.
В общем случае . Минимальное число узлов, необходимое для вычисления k-й производной, равно k + 1.
Оценка погрешности при численном дифференцировании может быть осуществлена по формуле
,
где n -- число узлов функции, k - порядок производной.
На практике чаще всего используются упрощенные формулы для равномерной сетки, при этом точность нередко повышается. Часто используются следующие формулы для трех узлов:
,
где h = x1 - x0 = const.
.
Исходя из общего вида интерполяционного полинома можно вывести формулы для более высокого порядка точности или для более высоких производных.
7 Численное решение систем линейных уравнений
Классы задач линейной алгебры
При численном решении большого круга задач в конечном итоге происходит их линеаризация, в связи с чем в соответствующих алгоритмах весьма широко используются методы линейной алгебры. В их числе:
· решение систем линейных алгебраических уравнений (СЛАУ);
· вычисление определителей матриц ;
· нахождение обратных матриц ;
· определение собственных значений и собственных векторов матриц ;
Постановка задачи решения СЛАУ:
, (1)
где -- квадратная матрица коэффициентов размерности n, - вектор неизвестных, -- вектор свободных коэффициентов. Иногда СЛАУ представляют в виде расширенной матрицы размерности n Ч n+1, где в качестве последнего столбца фигурирует вектор свободных коэффициентов-. В координатном представлении такая СЛАУ выглядит следующим образом:
. (2)
Для решения СЛАУ применяют в основном два класса методов: прямые (выполняемые за заранее известное количество действий) и итерационные (обеспечивающие постепенную сходимость к корню уравнения, зависящую от многих факторов). Прямые методы обычно применяются для решения систем порядка n < 200, для бульших n используются итерационные методы. Перед решением СЛАУ требуется проанализировать корректную постановку задачи:
1) Если - решение существует и единственно. Если же определитель равен нулю, то тогда, если матрица вырождена (т.е. ее можно преобразовать к виду, когда как минимум одна строка коэффициентов - нули) решений бесконечное множество, иначе решения не существует.
2) Если не имеет элементов с большими по модулю значениями - решение устойчиво (см. пример к главе 1). Показателем плохо обусловленных систем является .
Алгоритм метода Гаусса
1) Прямой ход.
Идея метода состоит в последовательном исключении неизвестных из системы n линейных уравнений. На примере первого уравнения системы (2) рассмотрим выражение для x1:
.
Подставим выражение для x1 во второе и все остальные уравнения системы:
.
Для расширенной матрицы коэффициентов это означает, что каждый элемент первой строки следует поделить на диагональный элемент, а все остальные строки преобразовать, как показано выше. Таким образом, станут равны нулю все коэффициенты первого столбца, лежащие ниже главной диагонали. Затем аналогичная процедура проводится со второй строкой матрицы и нижележащими строками, при этом первая строка и первый столбец уже не изменяются. И так далее до тех пор, пока все коэффициенты, лежащие ниже главной диагонали, не будут равны нулю.
Общие формулы прямого хода:
, (3)
k = 1…n, j = 1…n+1. Звездочкой отмечены элементы k-й строки с измененными значениями, которые будут подставлены в следующую формулу. Для определенности будем считать первый индекс - по строкам, второй - по столбцам.
, (4)
i = k +1…n, j = 1…n+1, k фиксировано в уравнении (3). Для уменьшения количества действий достаточно изменять значения элементов, находящихся выше главной диагонали.
2) Обратный ход.
Второй этап решения СЛАУ методом Гаусса называется обратным ходом и состоит в последовательном определении xk, начиная с xn, так как для последнего решение фактически получено. Общая формула:
. (5)
Таким образом, вычисление корней происходит за 2/3 n3 арифметических действий.
3) Выбор главного элемента.
Для уменьшения погрешности вычислений следует стремиться к тому, чтобы на главной диагонали матрицы стояли максимальные по модулю значения коэффициентов. Алгоритмически этого можно добиться, переставляя строки таким образом, чтобы на диагонали стоял наибольший по модулю элемент текущего столбца. Такая процедура называется выбором главного элемента и осуществляется всякий раз при переходе к новой строке в прямом цикле метода Гаусса.
4) Погрешность метода. Расчет невязок.
Точность результатов будет определяться только точностью выполнения арифметических операций при преобразовании элементов матрицы, т.е. ошибкой округления. Контроль правильности полученного решения осуществляется подстановкой полученных значений x1…xn в исходную систему уравнений и вычислением невязок, т.е. разностей между правыми и левыми частями уравнений:
, где k = 1…n. (6)
Специально отметим, что подставлять найденные значения следует в исходную (не преобразованную к верхнетреугольному виду) систему.
5) Преимущества и недостатки метода.
Преимущество метода в том, что он позволяет достичь результата за заранее известное и фиксированное число действий. Точность результатов будет определяться правильным выбором порядка коэффициентов в матрице и ее размерностью. Недостатком метода является резкое увеличение времени и погрешности вычислений с ростом n.
Блок-схема алгоритма метода Гаусса без выбора главного элемента.
Итерационные методы решения систем линейных уравнений.
Простейшим итерационным методом решения СЛАУ является метод простой итерации. При этом система уравнений (1) преобразуется к виду (2), а ее решение находится как предел последовательности (3), где {n} - номер итерации. Утверждается, что всякая система (2), эквивалентная (1), записывается в виде .
Теорема о достаточном условии сходимости метода простой итерации утверждает, что если норма матрицы (), то система уравнений (2) имеет единственное решение и итерационный процесс (3) сходится к решению со скоростью геометрической прогрессии.
Теорема о необходимом и достаточном условии сходимости метода простой итерации: Пусть система (2) имеет единственное решение. Итерационный процесс (3) сходится к решению системы (2) при любом начальном приближении тогда и только тогда, когда все собственные значения матрицы по модулю меньше 1.
На практике для обеспечения сходимости итерационных методов необходимо, чтобы значения диагональных элементов матрицы СЛАУ были преобладающими по абсолютной величине по сравнению с другими элементами.
Представим СЛАУ в следующей форме, удовлетворяющей (3):
(4)
Зададим начальные приближения и вычислим правую часть (4), получим новые приближения , которые опять подставим в систему (4). Таким образом организуется итерационный процесс, который обрывается по условию , где - заданная погрешность.
К ускорению сходимости приводит использование приближения к решениям путем последовательного уточнения компонентов, причем k-я неизвестная находится из k-го уравнения. Такая модификация итерационного метода носит название метода Зейделя:
Критерий сходимости метода Зейделя: Пусть -- вещественная симметричная положительно определенная матрица. Тогда метод Зейделя сходится.
Достоинствами метода простых итераций является простота программной реализации и более быстрый, по сравнению с линейными методами, поиск решения в матрицах большого размера. Недостатками являются сложный контроль условий сходимости и выбора начального приближения.
8)Интерполирование функций.
1) Необходимость: приблизить f(x) более простой функцией ф(х), совпадающей в узлах xi с f(xi), если f(x) определена только в узловых точках (результат эксперимента) или очень сложно вычисляется.
Условия Лагранжа : ф(х, с0, с1…сn) = fi,
0 <_i < n, где сi - свободные параметры, определяемые из данной системы уравнений.
С помощью интерполяции решают широкий круг задач численного анализа: дифференцирование и интегрирование функций, нахождение нулей и экстремумов, решение дифференцированных и т. д. Термин интерполяция употребляют, если х заключено между узлами, если он выходит за крайний узел, говорят об экстраполяции(при которой трудно гарантировать надежность приближения).
2) Пусть ф (х) = с0 + с1х + с2х2 +…+ сnxn (канонический вид полинома) ;сетка узлов может быть неравномерной.
Коэффициенты сi определяются из условий Лагранжа:
Получившаяся СЛАУ относительно свободных
параметров сi имеет решение, если среди узлов
хi нет совпадающих.Ее определитель - определитель Вандермонда:
Общая блок-схема:
3) Пусть задано n+1 значение функции f(x) в узлах xj
ф(х) = Pn(х) = i (x-xj)/(xi-xj) - полином Лагранжа.
Преимущества: потребуется решать СЛАУ для определения значения полинома в точке х.
Недостатки: для каждого х полином требуется читать заново.
Погрешность формулы: (*)
Увеличение числа узлов и, соответственно, степени полинома Pn(x) ведет к увеличению погрешности из-за роста производных .
4) ф(х) = Pn(x) = A0+A1(x-x0)+A2(x-x0)(x-x1)+…+An(x-x0)(x-x1)…(x-xn-1) - многочлен Ньютона для n+1 узла.
Коэффициенты Ф представляют собой разделенные разности и записываются в виде:
А0 = f0
A1 = (f0-f1)/(x0-x1) = f01
A2 = (f01-f02)/(x1-x2) = f012, где f02 = (f0-f2)/(x0-x2)
A3 = (f012-f013)/(x2-x3) = f0123 , где f013 = (f01-f03)/(x1-x3) , а f03 = (f0-f3)/(x0-x3)
и в общем случае Ak = (f01…k-1-f01…k)/(xk-1-xk)
Т.е. многочлен n-й степени выражается при помощи разделенных разностей через свои значения в узлах.
Преимущества: не решается СЛАУ, однако вычисление коэффициентов полинома не зависит от значения х и может быть вычислено только один раз. При добавлении нового узла также не происходит пересчета коэффициентов, кроме последнего.
После определения коэффициентов полинома Ньютона вычисление его значений при конкретных аргументах х наиболее экономично проводить по схеме Горнера:
P2(x) = A0+ (x-x0)(A1+(x-x2)(A3+…)…)
Погрешность определяется тем же соотношением (*)
Входящая в состав погрешности величина
(х-хi) = wn(x) ведет себя при постоянном шаге так, как показано на рисунке. Многочлен Ньютона имеет погрешность 0(hn+1) и обеспечивает n+1-й порядок точности интерполяции.
Между разделенными разностями и производными соответствующих порядков существует соотношение f <n>(x) ~ n! F01…n , где n - степень производной. Это используется в численном дифференцировании и при оценке погрешностей интерполяции.
Можно строить полиномы, не только проходящие через заданные точки, но и имеющие в них заданные касательные (интерполяционный многочлен Эрмита) или заданную кривизну. Количество всех полагаемых условий должно быть n-1, если n - степень полинома.
Основной недостаток интерполирования с помощью многочленов - неустранимые колебания, которые претерпевает кривая в промежутках между узлами.
При этом повышение степени интерполяционного полинома для большинства решаемых уравнений приводит не к уменьшению, а к увеличению погрешности.
Интерполяция сплайнами.
Происхождение термина “сплайны” связано с гибкой чертежной линейкой, которой пользовались для рисования гладких кривых, проходящих через заданные точки. Из теории упругости следует, что получающаяся кривая имеет постоянную кривизну и разрывы возникают лишь в третьей производной.
Обычно для сплайна выбирают кубический полином , определенный на интервале х из [xi-1, хi].
При этом вся кривая представляет собой набор таких кубических полиномов, с определенным образом подобранными коэффициентами аi , bi , ci , di , i- параметр сплайна.
Если вдоль сплайна совершается механическое движение, то непрерывность второй производной предполагает непрерывность ускорения и, следовательно, отсутствие резких изменений приложенной силы.
N+1 узлов
N интервалов
4N неизвестных
Условия подбора коэффициентов:
1)условия Лагранжа: ,
,
2)непрерывность первой и второй производной в узлах
фi'(xi) = фi+1'(xi); фi”(xi) = фi+1(xi)
3) условия в крайних узлах x0 и xn. Обычно задают условия свободных концов сплайна :
ф1”(x0) = 0, фn”(xn) = 0
Из полученных условий определяются зависимости между коэффициентами сплайнов:
В узле х = хi-1 коэффициент ai = fi-1.
В следующем узле x = xi выполняется условие
ai+bihi+cihi2+dihi3 = fi,
где элементарный шаг hi = xi - xi-1.
Потребуем непрерывности первой и второй производной на конце интервала
фi/(x) = bi+2ci(x-xi-1)+3di(x-xi-1)2 ,
фi//(x) = 2ci+6di(x-xi-1);
В узле x = xi первая производная
фi/(xi) = bi+2cihi+3dihi2 (1)
фi+1//(xi) = bi+1 (2)
Приравнивая (1) и (2), получаем bi +2cihi+3dihi2 = bi+1.
Вторая производная
фi//(xi) = 2ci+6cihi (3)
фi+1//(xi) = 2ci (4)
Приравнивая (3) и (4), получаем в свою очередь ci+3dihi = ci+1. Таким образом стыкуем все полиномы в узлах 1 ? i ? n-1. В крайних точках диапазона
ф1//(x0) = 2c1 = 0 > c1 = 0
ф1//(xn) = 2cn+6dnhn = 0 > cn +3dnhn = 0
Для всех 0 ? i ? n вышеприведенные соотношения представляют собой полную систему 4n линейных алгебраических уравнений относительно коэффициентов сплайнов, которую можно привести к системе ЛАУ, выразив коэффициенты ai, bi, di через ci и решить методом Гаусса или прогонки.
8 Система N линейных уравнений для коэффициентов сi:
для ,
где hi = xi-xi-1
После определения коэффициентов ci , 2N коэффициентов bi и di вычисляются по формулам:
,
И N уравнений для ,
Сплайновая интерполяция хороша тем, что требует знания в узлах только значений функции, но не ее производных.
Многомерная интерполяция
1) Последовательная интерполяция на прямоугольной сетке. Пусть заданы z i j = z(xi, yj) требуется найти z(x, y). Сначала при фиксированных yj0 найдем значение z(x, yj0),
Затем по полученному набору значений найдем z(x, y).
В случае интерполяции полиномом Лагранжа общая формула имеет вид
где k и m - количество узлов по сторонам прямоугольной сетки.
2) Треугольная конфигурация узлов.
z (x0, x1, y) = [z(x0, y)-z(x1, y)]/(x0-x1)
z (x, y0, y1) = [z(x, y0)-z(x,y1)]/(y0-y1)
Многочлен Лагранжева типа в этом случае имеет вид
9 Аппроксимация функций методом наименьших квадратов
Постановка задачи аппроксимации по МНК. Условия наилучшего приближения.
Если набор экспериментальных данных получен со значительной погрешностью, то интерполяция не только не требуется, но и нежелательна! Здесь требуется построить кривую, которая воспроизводила бы график исходной экспериментальной закономерности, т.е. была бы максимально близка к экспериментальным точкам, но в то же время была бы нечувствительна к случайным отклонениям измеряемой величины.
Введем непрерывную функцию ц(x) для аппроксимации дискретной зависимости f(xi), i = 0…n. Будем считать, что ц(x) построена по условию наилучшего квадратичного приближения, если
. (1)
Весу с для i-й точки придают смысл точности измерения данного значения: чем больше с, тем ближе аппроксимирующая кривая «притягивается» к данной точке. В дальнейшем будем по умолчанию полагать с = 1 для всех точек.
Рассмотрим случай линейной аппроксимации:
ц(x) = c0ц0(x) + c1ц1(x) + … + cmцm(x), (2)
где ц0…цm -- произвольные базисные функции, c0…cm - неизвестные коэффициенты, m < n. Если число коэффициентов аппроксимации взять равным числу узлов, то среднеквадратичная аппроксимация совпадет с интерполяцией Лагранжа, при этом, если не учитывать вычислительную погрешность, Q = 0.
Если известна экспериментальная (исходная) погрешность данных о, то выбор числа коэффициентов, то есть величины m, определяется условием:
. (3)
Иными словами, если , число коэффициентов аппроксимации недостаточно для правильного воспроизведения графика экспериментальной зависимости. Если , многие коэффициенты в (2) не будут иметь физического смысла.
Для решения задачи линейной аппроксимации в общем случае следует найти условия минимума суммы квадратов отклонений для (2). Задачу на поиск минимума можно свести к задаче поиска корня системы уравнений
, k = 0…m. (4).
Подстановка (2) в (1), а затем расчет (4) приведет в итоге к следующей системе линейных алгебраических уравнений:
Далее следует решить полученную СЛАУ относительно коэффициентов c0…cm. Для решения СЛАУ обычно составляется расширенная матрица коэффициентов, которую называют матрицей Грама, элементами которой являются скалярные произведения базисных функций и столбец свободных коэффициентов:
,
где , , j = 0…m, k = 0…m.
После того как с помощью, например, метода Гаусса найдены коэффициенты c0…cm, можно построить аппроксимирующую кривую или вычислить координаты заданной точки. Таким образом, задача аппроксимации решена.
Аппроксимация каноническим полиномом.
Выберем базисные функции в виде последовательности степеней аргумента x:
ц0(x) = x0 = 1; ц1(x) = x1 = x; цm(x) = xm, m < n.
Расширенная матрица Грама для степенного базиса будет выглядеть следующим образом:
.
Особенность вычислений такой матрицы (для уменьшения количества выполняемых действий) состоит в том, что необходимо сосчитать только элементы первой строки и двух последних столбцов: остальные элементы заполняются сдвигом предшествующей строки (за исключением двух последних столбцов) на одну позицию влево. В некоторых языках программирования, где отсутствует быстрая процедура возведения в степень, пригодится алгоритм расчета матрицы Грама, представленный далее.
Выбор базисных функций в виде степеней x не является оптимальным с точки зрения достижения наименьшей погрешности. Это является следствием неортогональности выбранных базисных функций. Свойство ортогональности заключается в том, что для каждого типа полинома существует отрезок [x0, xn], на котором обращаются в нуль скалярные произведения полиномов разного порядка:
, j ? k, с -- некоторая весовая функция.
Если бы базисные функции были ортогональны, то все недиагональные элементы матрицы Грама были бы близки к нулю, что увеличило бы точность вычислений, в противном случае при определитель матрицы Грама очень быстро стремится к нулю, т.е. система становится плохо обусловленной.
Аппроксимация ортогональными классическими полиномами.
Представленные ниже полиномы, относящиеся ко многочленам Якоби, обладают свойством ортогональности в изложенном выше смысле. То есть, для достижения высокой точности вычислений рекомендуется выбирать базисные функции для аппроксимации в виде этих полиномов.
1) Полиномы Чебышева.
Определены и ортогональны на [-1, 1] с весом . В интервал ортогональности всегда можно вписать область определения исходной функции с помощью линейных преобразований.
Строятся следующим образом (рекуррентная формула):
T0(x) = 1;
T1(x) = x;
Tk+1(x) = 2xTk(x) - Tk-1(x).
2) Полиномы Лежандра.
Определены и ортогональны на [-1, 1] с весом .
Строятся следующим образом (рекуррентная формула):
L0(x) = 1;
L1(x) = x;
.
Сглаживание и линейная регрессия.
Рассмотрим несколько наиболее простых с точки зрения программной реализации случаев аппроксимации (сглаживания).
1) Линейная регрессия.
В случае линейного варианта МНК (линейная регрессия) ц(x) = a + bx можно сразу получить значения коэффициентов a и b по следующим формулам:
,
,
где , .
2) Линейное сглаживание по трём точкам.
3) Линейное сглаживание по пяти точкам.
10 Решение нелинейных уравнений с одним неизвестным
Общие сведения о численном решении уравнений с одним неизвестным.
Пусть задана непрерывная функция f(x). Требуется найти корни уравнения f(x) = 0 численными методами - это и является постановкой задачи. Численное решение уравнения распадается на несколько подзадач:
1) Анализ количества, характера и расположения корней (обычно путем построения графика функции или исходя из физического смысла исследуемой модели). Здесь возможны следующие варианты:
· единственный корень;
· бесконечное множество решений;
· корней нет;
· имеется несколько решений, как действительных, так и мнимых (например, для полинома степени n). Корни четной кратности выявить сложно.
2) Локализация корней (разбиение на интервалы) и выбор начального приближения к каждому корню. В простейшем случае можно протабулировать функцию с заданным шагом. Если в двух соседних узлах функция будет иметь разные знаки, то между этими узлами лежит нечетное число корней уравнения (по меньшей мере один).
3) Вычисление каждого (или интересующего нас) корня уравнения с требуемой точностью. Уточнение происходит с помощью методов, изложенных ниже.
Метод дихотомии (бисекций).
Иначе называется методом половинного деления. Пусть задан начальный интервал [x0, x1], на котором f(x0)f(x1) ? 0 (т.е. внутри имеется не менее чем один корень). Найдем x2 = Ѕ (x0 + x1) и вычислим f(x2). Если f(x0)f(x2) ? 0, используем для дальнейшего деления отрезок [x0, x2], если > 0 - используем для дальнейшего деления отрезок [x1, x2], и продолжаем деление пополам.
Итерации продолжаются, пока длина отрезка не станет меньше 2о -- заданной точности. Тогда середина последнего отрезка даст значение корня с требуемой точностью. В качестве иного критерия можно взять | f(x)| ? оy.
Скорость сходимости метода невелика, однако он прост и надежен. Метод неприменим к корням четной кратности. Если на отрезке несколько корней, то заранее неизвестно, к какому из них сойдется процесс.
Если на заданном интервале предполагается несколько корней, то существует возможность последовательно исключать найденные корни из рассмотрения. Для этого воспользуемся вспомогательной функцией , где - только что найденный корень. Для функций f(x) и g(x) совпадают все корни, за исключением (в этой точке полюс функции g(x)). Для достижения требуемой точности рекомендуется грубо приблизиться к корню по функции g(x), а затем уточнить его, используя f(x).
Метод хорд.
Идея метода проиллюстрирована рисунком. Задается интервал [x0, x1], на котором f(x0)f(x1) ? 0, между точками x0 и x1 строится хорда, стягивающая f(x). Очередное приближение берется в точке x2, где хорда пересекает ось абсцисс. В качестве нового интервала для продолжения итерационного процесса выбирается тот, на концах которого функция имеет разные знаки. Условия выхода из итерационного цикла: или
| f(x)| ? оy.
Для вывода итерационной формулы процесса найдем точку пересечения хорды (описываемой уравнением прямой) с осью абсцисс: ax2 + b = 0, где ; b = f(x0) --- ax0.
Отсюда легко выразить .
Метод хорд в большинстве случаев работает быстрее, чем метод дихотомии. Недостатки метода те же, что и в предыдущем случае.
Метод Ньютона (касательных).
Пусть x0 - начальное приближение к корню, а f(x) имеет непрерывную производную. Следующее приближение к корню найдем в точке x1, где касательная к функции f(x), проведенная из точки (x0, f0), пересекает ось абсцисс. Затем точно так же обрабатываем точку (x1, f1), организуя итерационный процесс. Выход из итерационного процесса по условию .
Подобные документы
Традиционные языки высокоуровневого программирования. Обзор методов интегрирования. Оценка апостериорной погрешности. Численное решение систем линейных уравнений. Аппроксимация функций методом наименьших квадратов. Решение дифференциальных уравнений.
методичка [6,4 M], добавлен 23.09.2010Программирование нестандартных функций. Матрицы и операции над ними. Решение нелинейных уравнений и численное интегрирование. Оптимизация функции и численное дифференцирование. Аппроксимация функции: методы решения, описание программы, результаты.
курсовая работа [70,5 K], добавлен 12.08.2011Отделение корней методом простых интеграций. Дифференцирование и аппроксимация зависимостей методом наименьших квадратов. Решение нелинейного уравнения вида f(x)=0 методом Ньютона. Решение системы линейных уравнений методом Зейделя и методом итераций.
курсовая работа [990,8 K], добавлен 23.10.2011Матричная форма записи системы линейных уравнений, последовательность ее решения методом исключений Гаусса. Алгоритмы прямого хода и запоминания коэффициентов. Решение задачи о сглаживании экспериментальных данных с помощью метода наименьших квадратов.
курсовая работа [610,7 K], добавлен 25.06.2012Решение нелинейных уравнений методом простых итераций и аналитическим, простым и модифицированным методом Ньютона. Программы на языке программирования Паскаль и С для вычислений по вариантам в порядке указанных методов. Изменение параметров задачи.
лабораторная работа [191,0 K], добавлен 24.06.2008Особенности решения уравнений с одной переменной методом половинного деления. Оценка погрешности метода простой итерации. Суть решения уравнений в пакете Mathcad. Векторная запись нелинейных систем. Метод Ньютона решения систем нелинейных уравнений.
курсовая работа [2,1 M], добавлен 12.12.2013Решение уравнения методом половинного деления. Программа в Matlab для уравнения (x-2)cos(x)=1. Решение нелинейных уравнений методом Ньютона. Интерполяция заданной функции. Решение системы линейных алгебраических и обыкновенных дифференциальных уравнений.
курсовая работа [1,4 M], добавлен 15.08.2012Итерационные методы решения нелинейных уравнений, системы линейных алгебраических уравнений (СЛАУ). Решение нелинейных уравнений методом интерполирования. Программная реализация итерационных методов решения СЛАУ. Практическое применение метода Эйлера.
курсовая работа [1,6 M], добавлен 20.01.2010Решения алгебраических уравнений методом выделения корней. Аппроксимация функций методом наименьших квадратов; дихотомия, бисекция. Одномерная оптимизация многоэкстремальных функций; метод золотого сечения. Многомерная оптимизация градиентным методом.
курсовая работа [956,7 K], добавлен 04.03.2013Решение задачи Коши для дифференциального уравнения методом Рунге-Кутта и Адамса с автоматическим выбором шага и заданным шагом. Интерполирование табличной функции. Численное решение системы линейных алгебраических уравнений методами простой итерации.
методичка [35,8 K], добавлен 15.03.2009
