Курсовая работа Часть 1.Интерполяция Вариант 8. Вычислительная математика
Скачать 374.5 Kb.
|
Министерство образования и науки Российской Федерации Санкт-Петербургский государственный политехнический университет Институт Информационных технологий и управления Кафедра «Системный анализ и управление» Курсовая работа Часть 1.Интерполяция Вариант 8 дисциплина «Вычислительная математика» Выполнил: __________ Группа: __________ Проверил: __________ Санкт-Петербург 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. Методика экспериментов: 20 6.4. Теоретически ожидаемый результат. 20 6.5. Результаты исследования зависимости погрешности интерполирования от количества узлов интерполяции. 21 6.6. Результаты исследования сходимости. 22 ВведениеЗадача приближенного вычисления функции по известным значениям в заданных точках часто возникает в практической деятельности и теоретических задачах, например, при измерениях, определением экономических, медицинских, статистических показателей и др., когда нет математической модели исследуемого явления, а имеются лишь эмпирические данные, полученные в определенных точках, и возникает необходимость определения параметров в промежуточных точках. В этом случае прибегают к интерполяции. Кроме того, та же проблема имеет место, если заданная функция имеет сложный вид и ее вычисления сопряжены со значительными временными затратами. В этом случае при наличии удовлетворительной интерполяции выгоднее вычислить значение функции в некоторых узлах, а для промежуточных значений использовать интерполяцию. При этом выбирается функция таким образом, чтобы наилучшим образом удовлетворять имеющимся данным. Интерполяция применяется также в задаче обратного интерполирования, заключающаяся в нахождении обратной зависимости 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 = –1, c = –1, d = 1. 2. Анализ задания. Выбор способов решения задачи2.1. Аналитическое представление функцииНа рассматриваемом промежутке [c, d] интерполируемая функция кусочно-линейна, и на каждом линейном участке ее можно представить в виде y = kx + q, где k и q – постоянные. Видим два линейных участка: [c, 0] и [0, d]. На первом участке [–1, 0] функция не изменяется, значит k = 0, то есть y = q. При x = 0 y = 1, значит q =1, следовательно, на первом участке [–1, 0] уравнение функции имеет вид y = 1. На втором участке [0, 1] функция изменяется от 1 до –1, значит k = Δy/Δx = =(–1–1)/(1–0)= –2, то есть y = –2x + q. При x = 0 y = 1, значит q = 1, итак, на втором участке [0, 1] уравнение функции имеет вид y= –2x + 1. Таким образом, интерполируемая функция имеет вид: (1) 2.2. Исходные данныеВ соответствии с заданием, необходимо сформировать таблицу с исходными данными для равномерных разбиений и разбиений по Чебышеву. Для этого на языке программирования MatLab была подготовлена функция Slope(x), рассчитываемая по формуле (1). Код функции приведен в разделе 4.2 в составе листинга программы. Расчет выполнялся в следующем порядке: 1) заданы значения c и d: с= –1; d=1; 2) задано значение n, например, n=2; 3) определен массив i: i=1:n+1; 4) определен массив x: а) равномерное разбиение: шаг h=(d–c)/n, xi =c+(d–c)(i–1)/n: 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=slope(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. Формирование массива интерполируемой функцииИнтерполируемая функция была реализована в виде подпрограммы Slope(x) (см. раздел 4.2). Алгоритм ее вычисления следующий: n := length(x) дляk = 1..n Массив интерполируемой функции Fn реализован следующим образом: x1=c:(d-c)/m:d; Fn=Slope(x1); 3.3. Формирование массива значений интерполяционного полиномаОпределение узлов описано в разделе 2.2. Алгоритм формирования массива значений интерполяционного полинома приведен ниже: д ля xi = c + (d – c) * (i – 1) / m yk = Newton(n, xk) Здесь Newton(n, x) – подпрограмма вычисления интерполяционного полинома n-го порядка. Он рассчитывается по формулам (3–4) по следующему алгоритму: Формируется матрица разделенных разностей fij по формуле (4). Создается заготовка матрицы: f=zeros(n+1, n+1). Заполняется первыйстолбец (значения функции в точках xk): f(:, 1)=y’. Остальные столбцы формируются в цикле: д ля д ля ; F=f(1,:); Построение полинома: д ля x1(i)=c+i*(d–c)/m задается начальное значение при формировании : – первое значение в сумме д ля –формирование сомножителя (слагаемого) – формирование суммы полинома 3.4. Построение графиковГрафики строятся по полученным значениям Fn, Newton и y. 4. Разработка программы4.1. Структура программыВ соответствии с разработанным алгоритмом, программа состоит из следующих подпрограмм: Интерполируемая функция Slope(x). Основной модуль Kurs. Модуль состоит из следующих блоков: - инициализация; - расчет разбиения и значений в узлах; - расчет матрицы разделенных разностей; - расчет полинома; - расчет исходной функции; - построение графика. 4.2. Листинг программыfunction Kurs() n=40; % число узлов c=-1; % начало 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=Slope(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; Newton(i)=F(1); for j=2:n+1 w=w*(xi(i)-x(j-1)); Newton(i)=Newton(i)+w*F(j); end end Fn=Slope(xi); % Вызов функции subplot(2,1,1) plot(x,y,'s','MarkerFaceColor','c','MarkerSize',8); hold on; grid on; plot(xi,Newton,'b','LineWidth',3); plot(xi,Fn,'g','LineWidth',2); title 'Узлы, исходная функция и полином' subplot(2,1,2) hold on; plot(xi,(Fn-Newton),'r'), grid on; axis tight; title 'Погрешность' for i=1:m if xi(i) D1(i)=abs(Fn(i)-Newton(i)); elseif xi(i) D2(i)=abs(Fn(i)-Newton(i)); else D3(i)=abs(Fn(i)-Newton(i)); end D(i)=abs(Fn(i)-Newton(i)); end [max(D); max(D1); max(D2); max(D3)] %===< Интерполируемая функция >===% function [Y]=Slope(x) for i=1:length(x) if x(i)<-1 Y(i)=0; elseif x(i)<0 Y(i)=1; elseif x(i)<=1 Y(i)=-2*x(i)+1; else Y(i)=0; end end 5. ТестированиеЦелью тестирования является подтверждение того, что написанная программа работает в соответствии с техническим заданием. Для того, чтобы доказать, что программа функционирует так, как положено, необходимо, чтобы соблюдалось выполнение тестового критерия. Используем в качестве тестового критерия следующее: Значения исходной функции и полученного интерполирующего полинома должны совпадать в узлах интерполирования при любом количестве узлов. Для проверки выполнения тестового критерия рассмотрим несколько случаев с различным числом узлов: n=3, 7, 15, 41 (степени многочлена 2, 6, 14, 40). Графические результаты тестирования приведены на рисунках 2–9. 1. Равномерное разбиение.
2. Разбиение по Чебышеву.
Из графиков видно, что интерполируемая функция и ее интерполяция совпадают в узлах интерполирования. Это означает, что выполнен тестовый критерий, что подтверждает правильность подготовленной программы. 6. Исследования6.1. Цель исследования:1. Исследование зависимости погрешности интерполирования от количества узлов. 2. Исследование сходимости метода интерполяции. 6.2. Программа исследования:1. Определение зависимости погрешности интерполирования от количества узлов интерполяции путем определения наибольшего отклонения интерполирующего полинома от исходной функции для равномерного разбиения и разбиения по Чебышеву для степеней многочлена n = 2, 6, 14, 40. 2. Определение числа точек разбиения, при котором начинается расхождение, либо принятие решения об отсутствии расхождения интерполяции. 6.3. Методика экспериментов:Необходимые параметры эксперимента задаются в функции Kurs путем редактирования. Это число узлов n (строка 3), число точек, используемых для построения векторов и графиков m (строка 6) и определение стратегии выбора узлов (строки 9, 10). Последнее состоит в задании разбиения x(i) путем задействования одной строки, задающей разбиение (по Чебышеву или равномерное разбиение) и закрытие комментарием другой строки. После задания параметров эксперимента в функции Kurs она запускается из командного окна MatLab. При этом выполняется расчет и вывод наибольшего отклонения интерполируемой функции от интерполирующего многочлена (т.е. абсолютной погрешности) в виде: max(abs(Fn-Newton)) %<--- Вывод максимальной погрешности Здесь массив Fn содержит значения заданной функции Slope(x), а массив Newton – соответствующие им рассчитанные значения интерполирующего полинома. Последовательно устанавливая n и стратегию выбора узлов в соответствии с программой исследования, получаем графики и наибольшую для заданных параметров величину погрешности. Таким образом, выполняется исследование погрешности от n и стратегии выбора узлов. Исследования сходимости выполняется аналогично, только величина n последовательно увеличивается до 40 и определяется, возникает ли рост погрешности, и если возникает, то при каком n. 6.4. Теоретически ожидаемый результат.Теоретически выбор в качестве узлов интерполирования разбиения по Чебышеву является наилучшим из данных двух. При интерполяции с равноотстоящими узлами ошибка интерполяции должна быть существенно выше на краях промежутка, а при интерполяции с чебышевскими узлами ошибка более равномерно распределена на интервале интерполирования. Как было отмечено во введении, для любой функции с ограниченными производными на заданном промежутке при выборе для интерполяции чебышевских узлов метод интерполяции сходится. Заданная функция, однако, имея излом в нуле, содержит в производных, начиная со второй, дельта-функцию, следовательно, ее производные не ограничены, и условие теоремы о сходимости интерполяции при разбиении по Чебышеву не выполняется. 6.5. Результаты исследования зависимости погрешности интерполирования от количества узлов интерполяции.Полученные значения погрешности (наибольшего различия интерполируемой функции и интерполирующего полинома) сведены в таблицу:
По результатам тестирования можно сделать вывод о том, что при любом выборе числа узлов метод Чебышева дает более точный результат по сравнению с методом равномерных узлов (погрешность при всех значениях n меньше). По графикам (раздел 5) видно, что при интерполяции с равноотстоящими узлами ошибка интерполяции значительно выше на краях промежутка интерполяции:
При интерполяции с чебышевскими узлами ошибка более равномерно распределена на интервале интерполирования:
Для всех методов первоначальное увеличение числа точек разбиения ведет к увеличению точности результата, однако при последующем увеличении этого числа точность падает, причем возникновение этой расходимости при равномерном выборе точек разбиения возникает раньше (рис. 2–9). При этом интерполирующий полином, как правило, имеет экстремумы, точки перегиба на линейных участках интерполируемой функции, но не имеет излома в точке излома исходной функции, то есть производные у него всюду непрерывны. 6.6. Результаты исследования сходимости.Сходимость интерполяции означает неограниченное приближение интерполирующего полинома к интерполируемой функции при неограниченном увеличении числа узлов интерполяции. Строгое доказательство сходимости может быть получено только аналитически. Численный эксперимент может установить только расходимость. При этом причина расходимости (погрешность метода или погрешность вычислений) опять же может быть установлена только при помощи дополнительного анализа. В связи с этим в целях исследования сходимости выполнялось определение наибольшей погрешности (разницы) построенного полинома от исследуемой функции при различной стратегии выбора узлов по тестовому равномерному разбиению с количеством точек m=100. Для этого аналогично предыдущему разделу задавалась степень полинома (начиная от 2) и стратегия выбора узлов (по Чебышеву или равномерно) и осуществлялся расчет максимальной погрешности на исследуемом промежутке. Исследование показало, что увеличения числа точек разбиения приводит к расходимости для равномерного разбиения, начиная уже с n=6. Для разбиения по Чебышеву эксперимент проводился до n=40, и расхождения в этих пределах обнаружено не было. В то же время такой метод исследования не может гарантировать, что при дальнейшем увеличении степени многочлена не будет возникать расходимости, связанной, например, с неточностью вычислений. ЗаключениеВ ходе выполнения работы изучены основные методы интерполяции. Разработан алгоритм интерполяции заданной кусочно-линейной функции по формуле Ньютона и подготовлена программа в математической системе MatLab, реализующая разработанный алгоритм. Исследована зависимость погрешности интерполяции от степени интерполирующего полинома, а также исследована сходимость интерполяции. Установлено, что погрешность монотонно убывает для равномерного разбиения – до степени полинома n = 5 включительно, для разбиения по Чебышеву вплоть до n = 40 расходимости не обнаружено. Список используемой литературыИльина В.А., Силаев П.К. Численные методы для физиков-теоретиков. I. — Москва-Ижевск: Институт компьютерных исследований, 2003. – 132 с. Лапчик М.П. Численные методы: Учеб. пособие для студ. вузов. – М., 2004 – 384 с. Калиткин Н.Н. Численные методы. – М.: Наука, 1978. – 512 с. Кирсяев А.Н., Куприянов В.Е. Вычислительная математика. Конспект лекций. – СПб., 2011. – 173 с. Турчак Л. И., Плотников П. В. Основы численных методов: Учебное пособие. — 2-е изд., перераб. и доп. — М.: ФИЗМАТЛИТ, 2003. — 304 с. |