Курсовая работа Часть 1.Интерполяция Вариант 27. Курсовая. Вычислительная математика
Скачать 383.5 Kb.
|
Министерство образования и науки Российской Федерации Санкт-Петербургский государственный политехнический университет Институт Информационных технологий и управления Кафедра «Системный анализ и управление» Курсовая работа Часть 1.Интерполяция Вариант 27 дисциплина «Вычислительная математика» Выполнил: __________ Группа: __________ Проверил: __________ Санкт-Петербург 2016 Содержание 2.1. Аналитическое представление функции 8 2.2. Исходные данные 8 2.3. Выбор метода 10 2.4. Расчетные формулы метода 10 3.1. Исходные данные 10 3.2. Формирование массива интерполируемой функции 11 3.3. Формирование массива значений интерполяционного полинома 11 3.4. Построение графиков 12 4.1. Структура программы 12 4.2. Листинг программы 12 6.1. Цель исследования: 19 6.2. Программа исследования: 19 6.3. Методика экспериментов: 19 6.4. Теоретически ожидаемый результат. 20 6.5. Результаты исследования зависимости погрешности интерполирования от степени интерполирующего полинома. 20 6.6. Результаты исследования сходимости. 21 ВведениеПроблема приближенного вычисления функции по известным значениям в заданных точках часто возникает в задачах, связанных с измерениями, определением экономических, медицинских, статистических показателей и др., когда нет математической модели исследуемого явления, а имеются лишь эмпирические данные, полученные в некоторых точках, и возникает необходимость определения параметров в промежуточных точках. В этом случае прибегают к интерполяции. Кроме того, та же проблема имеет место, если заданная функция имеет сложный вид и ее вычисления сопряжены со значительными временными затратами. В этом случае при наличии удовлетворительной интерполяции выгоднее вычислить значение функции в некоторых узлах, а для промежуточных значений использовать интерполяцию. При этом выбирается функция таким образом, чтобы наилучшим образом удовлетворять имеющимся данным. Интерполяция применяется также в задаче обратного интерполирования, заключающаяся в нахождении обратной зависимости x(y) по известным точкам y(xi). Примером обратного интерполирования может служить задача о нахождении корней уравнения. Перечисленные примеры задач, часто встречающихся на практике, показывают актуальность проблемы интерполяции функции по известным значениям. Таким образом, задача интерполяции функции f(x) на интервале [a, b] состоит в нахождении коэффициентов ci так, чтобы отклонение заданной функции f(x) от обобщенного полинома (1) было минимальным на интервале [a, b]. Здесь φi(x) − набор действительных линейно независимых (чебышевских) функций, при помощи которых строится интерполяция. Решение такой задачи сводится к решению системы линейных алгебраических уравнений (СЛАУ). Действительно, из (1) для набора точек xk, k = 0..n, получим систему представляющей собой СЛАУ относительно ci. Для ее однозначного решения необходимо, чтобы определитель системы (2) был отличен от нуля (при этом предполагается, что все xk не равны между собой, т.е. функция задана в различных точках). В качестве функций φi(x), используемых для аппроксимации, как правило, выбирают: 1) полиномиальные функции 1, x, x2, … либо их линейные комбинации (полином Лагранжа, полином Ньютона) − алгебраическое интерполирование; 2) тригонометрические функции 1, sin(x), cos(x), sin(2x), cos(2x), … − тригонометрическое интерполирование; 3) exp(a0x), exp(a1x), exp(a2x), … либо их линейные комбинации − экспоненциальное интерполирование; 4) различные специальные функции − полиномы Лежандра, Эрмита, Чебышева, Ляггера и др. Отдельное место занимает рациональная интерполяция, при которой строится дробно-рациональная функция, представляющая собой отношение , равное f(x) в заданных точках. Основными проблемами численных методов при интерполяции являются погрешность и сходимость. Исследуя сходимость метода интерполяции, будем говорить, что при заданной стратегии выбора узлов метод интерполяции сходится, если ΔPn= при n . При n→∞ в случае использования чебышевской сетки ΔPn→0; для равномерной сетки в аналогичных условиях ΔPn→∞. Рассмотренная выше ситуация обобщается следующими теоремами: Теорема 1. Какова бы не была стратегия выбора узлов, найдется непрерывная на интервале [a, b] функция такая, что метод интерполяции расходится, то есть En(y)→∞ при n→∞. Однако для функций, обладающих ограниченными производными, стратегия, обеспечивающая сходимость, существует. Теорема 2. Пусть в качестве узлов интерполяции выбирают чебышевские узлы. Тогда для любой функции, обладающей ограниченными производными на [a, b], метод интерполяции сходится. Погрешность интерполирования ε = f(x)−φ(x) зависит: от исходных данных, в частности, от числа узлов; от расположения узлов xi, от выбранного правила интерполирования (выбора класса функций φ). На практике имеются два пути приближения к пределу. Первый состоит в уменьшении шага сетки, второй − в увеличении степени полинома. При этом если первый способ всегда приводит к увеличению точности для почти любой интерполируемой функции, то увеличение степени многочлена приводит к улучшению результата далеко не всегда. Например, степенной многочлен можно рассматривать как ряд Тейлора, то есть коэффициенты при степенях x есть производные высших степеней, но не все функции имеют отличные от нуля производные высших порядков. В этом случае введение высших степеней бессмысленно. Поэтому при большом количестве заданных точек используют сплайны − участки интерполяции, строящиеся на нескольких ближайших точках (на практике, от 2 до 4). Оценивать ошибку интерполяции совершенно необходимо, иначе можно получить бессмысленный результат. Ошибки при полиномиальной интерполяции, как правило, возрастают вблизи концов интервала, на котором лежат x1..xn. Кроме того, обычно с возрастанием количества точек n точность растет. Однако это верно только для тех функций, которые хорошо аппроксимируются полиномом, т.е. для функций с большим радиусом сходимости степенного ряда. То есть, если полюса функции лежат достаточно далеко от точек x1..xn, то интерполяция тем точнее, чем больше количество точек n. 1. Постановка задачиСравнить графики заданной функции f(x) и интерполяционных полиномов Pn(x) для n=2,6,14, на интервале [c, d] для двух вариантов выбора узлов: а) равномерно, с шагом h=(d−c)/n; б) по Чебышеву. Число узлов n+1, n − степень интерполяционного полинома. Для каждого способа выбора узлов рассчитать погрешность интерполирования в зависимости от числа узлов, сделать заключение о сходимости метода интерполяции. Исходные данные задачи: Интерполируемая функция задана графиком (рис. 1) с параметрами: a=1, b= –2, c= –2, d=1.
2. Анализ задания. Выбор способов решения задачи2.1. Аналитическое представление функцииНа рассматриваемом промежутке [c, d] функция является кусочно-линейной, поэтому для каждого линейного участка ее можно представить в виде y=kx+p, где k и p – постоянные. Видим четыре линейных участка: [–2, –1], [–1, 0], [0, 1] и [1, 2]. На первом участке [–2, –1] функция изменяется от 0 до 1, значит k=Δy/Δx= =(1-0)/(–1–(–2))=1, то есть y=x+p. При x= –2 y=0, значит p=2, итак, на первом участке [–2, –1] уравнение функции: y=x+2. На втором участке [–1, 0] функция изменяется от 1 до 0, значит k=Δy/Δx= =(1–0)/(0–(–1))=–1, то есть y=–x+p. При x=0 y=0, значит p=0, итак, на втором участке [–1, 0] уравнение функции y= –x. На третьем участке [0, 0.5] функция изменяется от 0 до –2, значит k=Δy/Δx= =(–2–0)/(0.5–0)= –4, то есть y= –4x+p. При x=0 y=0, значит p=0, итак, на третьем участке [0, 0.5] уравнение функции: y= –4x. На четвертом участке [0.5, 1] функция изменяется от –2 до 0, значит k=Δy/Δx=(0–(–2))/(1–0.5)=4, то есть y= 4x+p. При x=1 y=0, значит p= –4, итак, на четвертом участке [0.5, 1] уравнение функции y=4x–4. Таким образом, в аналитическом виде интерполируемая функция запишется в следующей форме: (1) 2.2. Исходные данныеСформируем таблицу с исходными данными для равномерных разбиений и разбиений по Чебышеву. Для этого на языке MatLab подготовим функцию updown(x) в соответствии с формулой (1). Код функции приведен в листинге программы в разделе 4.2. Расчет данных выполнялся следующим образом: 1) были заданы значения c и d: с= –2; d=1; 2) задавалось значение n, например, n=2; 3) задавался массив i: i=1:n+1; 4) задавался массив x: а) равномерное разбиение: шаг h=(d–c)/n, xi=c+(i – 1)*h, x=c:(d–c)/n:d б) разбиение по Чебышеву: : (2) x=(c+d)/2 + (d–c)/2 * cos((2*i–1)./(2*(n+1))*pi) 5) определялся массив yi: y=updown(x) Вычисленные данные для n=2, 6, 14 внесены в таблицы. Равномерное разбиение. n=2:
n=6:
n=14
разбиение по Чебышеву. n=2:
n=6:
n=14
2.3. Выбор методаНеобходимо выполнить алгебраическую интерполяцию функции в виде полинома Лагранжа или Ньютона. Для расчетов на ЭВМ удобнее использовать формулу Ньютона, поэтому остановим свой выбор на этом методе. 2.4. Расчетные формулы методаИнтерполяционным многочленом Ньютона называется многочлен (3) где , а разделенные разности f(x0, x1, …, xk) могут быть определены по рекуррентной формуле , (4) , ; . При этом для всех i, j должно выполняться соотношение xj– xi 0, что, справедливо как для разбиения Чебышева, так и для равномерного разбиения. 3. Разработка алгоритма3.1. Исходные данныеИсходными данными в данной задаче являются: постоянные a, b, c, d; степень интерполяционного полинома n; количество точек на графике m. 3.2. Формирование массива интерполируемой функцииИнтерполируемая функция оформлена в виде подпрограммы с именем updown(x). Ее листинг приведен в разделе 4.2. Массив интерполируемой функции UD формировался следующим образом: x1=c:(d-c)/m:d; UD=updown(x1); Алгоритм вычисления следующий (предполагается, что аргумент x может быть вектором): N:=length(x) дляk=1..N 3.3. Формирование массива значений интерполяционного полиномаПорядок расчета узлов описан в разделе 2.2. Алгоритм формирования массива значений интерполяционного полинома следующий: д ля xi=c+(d–c)*(i–1)/m yk=P (xk) Здесь Pn(x) – подпрограмма вычисления интерполяционного полинома n-го порядка. Его расчет производится по формулам (3, 4). Алгоритм расчета: Формирование матрицы fij разделенных разностей по формуле (4). Заготовка матрицы, заполненная нулями: f=zeros(n+1, n+1). Первый столбец (интерполируемая функция в точках xk): f(:, 1)=y’. Цикл формирования оставшихся столбцов: д ля д ля ; F=f(1,:); Построение полинома: д ля xi(i) = c+i*(d–c)/m начальное значение при формировании : – первое значение в сумме Д ля –формирование сомножителя (слагаемого) – формирование суммы полинома 3.4. Построение графиковГрафики строятся по вычисленным векторам UD, P и y. 4. Разработка программы4.1. Структура программыВ соответствии с разработанным алгоритмом, программа состоит из следующих подпрограмм: Интерполируемая функция updown(x). Основной модуль curs. Он состоит из следующих блоков: - блок инициализация; - блок расчета разбиения - блок расчета значений в узлах разбиения; - блок расчета матрицы разделенных разностей; - блок расчета значений интерполирующего полинома; - блок расчета значений заданной функции; - блок построения графиков; - блок вычисления максимальных погрешностей. 4.2. Листинг программыfunction Curs() n = 50; %// число узлов c = -2; %// начало d = 1; %// конец m = 100; %// число точек %// i = 1 : n + 1; %// номера узлов x(i) = (c + d)/2 + (d - c)/2 * cos(pi * (2*i - 1)/2/(n + 1)); %// по Чебышеву %x(i) = c + (d - c)*(i - 1)/n; %// Равномерное y=updown(x); %// вызов функции %// f = zeros(n + 1, n + 1); f(:, 1) = y'; % 1-й столбец for i = 2 : n + 1 for j = 1 : n - i + 2 f(j, i) = (f(j + 1, i - 1) - f(j, i - 1))/(x(j + i - 1) - x(j)); end; end; F = f(1, :); xi = c : (d - c)/m : d; for i = 1 : m + 1 w = 1; P(i) = F(1); for j = 2 : n + 1 w = w*(xi(i) - x(j - 1)); P(i) = P(i) + w*F(j); end end UD = updown(xi); % Вызов функции subplot(2, 1, 1) plot(x, y, 's', 'MarkerFaceColor', 'c', 'MarkerSize', 8); hold on; grid on; plot(xi, P, 'b', 'LineWidth', 3); plot(xi, UD, 'g', 'LineWidth', 2); title 'Исходная функция, полиномная интерполяция и узлы' subplot(2, 1, 2) hold on; plot(xi, (UD - P), 'r'), grid on; axis tight; title 'Погрешность' for i = 1 : m if xi(i) < -1 R1(i) = abs(UD(i) - P(i)); elseif xi(i) < 0 R2(i) = abs(UD(i) - P(i)); else R3(i) = abs(UD(i) - P(i)); end R(i) = abs(UD(i) - P(i)); end [n; max(R); max(R1); max(R2); max(R3)] %//***< Интерполируемая функция >***//% function [y]=updown(x) N = length(x); for k=1:N if (x(k) < -2) y(k) = 0; elseif x(k) <= -1 y(k) = x(k) + 2; elseif x(k) <= 0 y(k) = -x(k); elseif x(k) <= 0.5 y(k) = -4 * x(k); elseif x(k) <= 1 y(k) = 4 * x(k) - 4; else y(k) = 0; end end 5. ТестированиеТестирование выполняется для подтверждение того, что написанная программа работает в соответствии с техническим заданием. Для доказательства этого необходимо выполнение тестового критерия. Используем в качестве тестового критерия следующее: Значения исходной функции и полученного интерполирующего полинома должны совпадать в узлах интерполирования при любом количестве узлов. Для проверки выполнения тестового критерия рассмотрим несколько случаев с различным числом узлов: n = 3, 7, 15, 51 (степени многочлена 2, 6, 14, 50). Графические результаты тестирования приведены на рисунках 2–9. 1. Равномерное разбиение.
2. Разбиение по Чебышеву.
Из графиков видно, что интерполируемая функция и ее интерполяция совпадают в узлах интерполирования. Это означает, что выполнен тестовый критерий, что подтверждает правильность подготовленной программы. 6. Исследования6.1. Цель исследования:1. Исследование зависимости погрешности интерполирования от количества узлов. 2. Исследование сходимости метода интерполяции. 6.2. Программа исследования:1. Определение зависимости погрешности интерполирования от количества узлов интерполяции путем определения наибольшего отклонения интерполирующего полинома от исходной функции для равномерного разбиения и разбиения по Чебышеву для степеней многочлена n=2, 6, 14, 50. 2. Определение числа точек разбиения, при котором начинается расхождение, либо принятие решения об отсутствии расхождения интерполяции. 6.3. Методика экспериментов:Необходимые параметры эксперимента задаются в функции Interpolation путем редактирования. Это число узлов n (строка 3), число точек, используемых для построения векторов и графиков m (строка 4) и определение стратегия выбора узлов (строки 9, 10). Последний состоит в задании разбиения x(i) путем задействования одной строки, задающей разбиение (по Чебышеву или равномерное разбиение) и закрытие комментарием другой строки. После задания параметров эксперимента в функции Interpolation она запускается из командного окна MatLab. При этом выполняется расчет и вывод наибольшего отклонения интерполируемой функции от интерполирующего многочлена (т.е. абсолютной погрешности) в виде: max(abs(UD-P)) %<--- Вывод максимальной погрешности Здесь массив UD содержит значения заданной функции updown(x), а массив P – соответствующие им рассчитанные значения интерполирующего полинома. Последовательно устанавливая n и стратегию выбора узлов в соответствии с программой исследования, получаем графики и наибольшую для заданных параметров величину погрешности. Таким образом, выполняется исследование погрешности от n и стратегии выбора узлов. Исследования сходимости выполняется аналогично, только величина n последовательно увеличивается до 50 и определяется, возникает ли рост погрешности, и если возникает, то при каком n. 6.4. Теоретически ожидаемый результат.Теоретически выбор в качестве узлов интерполирования разбиения по Чебышеву является наилучшим из данных двух. При интерполяции с равноотстоящими узлами ошибка интерполяции должна быть существенно выше на краях промежутка, а при интерполяции с чебышевскими узлами ошибка более равномерно распределена на интервале интерполирования. Как было отмечено во введении, для любой функции с ограниченными производными на заданном промежутке при выборе для интерполяции чебышевских узлов метод интерполяции сходится. Заданная функция, однако, имея излом в нуле, содержит в производных, начиная со второй, дельта-функцию, следовательно, ее производные не ограничены, и условие теоремы о сходимости интерполяции при разбиении по Чебышеву не выполняется. 6.5. Результаты исследования зависимости погрешности интерполирования от степени интерполирующего полинома.Полученные значения погрешности (наибольшего различия интерполируемой функции и интерполирующего полинома) сведены в таблицу:
По результатам тестирования можно сделать вывод о том, что при любом выборе числа узлов метод Чебышева дает более точный результат по сравнению с методом равномерных узлов (погрешность при всех значениях n меньше). По графикам (раздел 5) видно, что при интерполяции с равноотстоящими узлами ошибка интерполяции значительно выше на краях промежутка интерполяции (за исключением полинома 2 степени):
При интерполяции с чебышевскими узлами ошибка более равномерно распределена на интервале интерполирования:
Для всех методов первоначальное увеличение числа точек разбиения ведет к увеличению точности результата, однако при последующем увеличении этого числа точность падает, причем возникновение этой расходимости при равномерном выборе точек разбиения возникает раньше (рис. 2–9). При этом интерполирующий полином, как правило, имеет экстремумы, точки перегиба на линейных участках интерполируемой функции, но не имеет излома в точке излома исходной функции, то есть производные у него всюду непрерывны. 6.6. Результаты исследования сходимости.Сходимость интерполяции означает неограниченное приближение интерполирующего полинома к интерполируемой функции при неограниченном увеличении числа узлов интерполяции. Строгое доказательство сходимости может быть получено только аналитически. Численный эксперимент может установить только расходимость. При этом причина расходимости (погрешность метода или погрешность вычислений) опять же может быть установлена только при помощи дополнительного анализа. В связи с этим в целях исследования сходимости выполнялось определение наибольшей погрешности (разницы) построенного полинома от исследуемой функции при различной стратегии выбора узлов по тестовому равномерному разбиению с количеством точек m=100. Для этого аналогично предыдущему разделу задавалась степень полинома (начиная от 2) и стратегия выбора узлов (по Чебышеву или равномерно) и осуществлялся расчет максимальной погрешности на исследуемом промежутке. Исследование показало, что увеличения числа точек разбиения приводит к расходимости для обоих стратегий выбора узлов, но для равномерного разбиения расходимость начинается значительно раньше. Монотонное увеличение точности происходит: - для равномерного разбиения – до n = 4 включительно (погрешность 0.7174); - для разбиения по Чебышеву – до n = 9 включительно (погрешность 0.2528). Далее при увеличении n имеется участок неустойчивости, когда увеличение числа узлов приводит то к увеличению точности, то к ее уменьшению, после чего начинается резкий рост погрешности: - для равномерного разбиения – при n=26; - для разбиения по Чебышеву – при n=45. Так как теоретически разбиение по методу Чебышева должно приводить к сходимости, полученная расходимость может быть объяснена неточностью вычислений. ЗаключениеВ ходе выполнения работы были изучены основные методы интерполяции. Разработан алгоритм интерполяции кусочно-линейной функции по формуле Ньютона, а также программа в математической системе MatLab, которая реализует этот алгоритм. Исследована зависимость погрешности интерполяции от степени интерполирующего полинома, а также исследована сходимость степенной интерполяции. Получен результат, что погрешность монотонно убывает для равномерного разбиения – до степени полинома n=5 включительно (погрешность 0.3247), для разбиения по Чебышеву – до n=8 включительно (погрешность 0.2196). Далее имеется участок неустойчивости, когда увеличение числа точек разбиения приводит то к увеличению, то к уменьшению точности, после чего начинается резкий рост погрешности: для равномерного разбиения – при n=26; для разбиения по Чебышеву – при n=45. Полученная расходимость интерполяции при использовании разбиения по Чебышеву, по всей вероятности, связана с неточностью вычислений. Список используемой литературыИльина В.А., Силаев П.К. Численные методы для физиков-теоретиков. I. — Москва-Ижевск: Институт компьютерных исследований, 2003. – 132 с. Лапчик М.П. Численные методы: Учеб. пособие для студ. вузов. – М., 2004 – 384 с. Калиткин Н.Н. Численные методы. – М.: Наука, 1978. – 512 с. Кирсяев А.Н., Куприянов В.Е. Вычислительная математика. Конспект лекций. – СПб., 2011. – 173 с. Иванов А.П., Олемской И.В. Численные методы. Часть II / учеб. пособие. – СПб., 2012. – 80 с. |