Презентация по дисциплине Компьютерное моделирование Аппроксимация и интерполяция функций
Аппроксимация, интерполяция, экстраполяция Аппроксимация (приближение) функции f(x) - нахождение такой функции g(x) (аппроксимирующей функции), которая была бы близка к заданной. Интерполяция – определение промежуточных значений функции fj по известному набору значений функции fi. Экстраполяция – определение значений функции за пределами первоначально известного интервала. f1 x1 f2 x2 … … fi xi … … fn xn
Интерполяция Интерполяция имеет как практическое, так и теоретическое значение. На практике часто возникает задача о восстановлении непрерывной функции по ее табличным значениям, например полученным в ходе некоторого эксперимента. Для вычисления многих функций оказывается эффективно приблизить их полиномами или дробно-рациональными функциями. Теория интерполирования используется для численного интегрирования, для получения методов решения дифференциальных и интегральных уравнений.
Задача интерполяции Пусть функция f(x) задана таблицей своих значений xi, yi: на интервале [a; b]: Задача интерполяции - найти функцию g(x), принимающую в точках xi те же значения yi точки xi – узлы интерполяции условие g(x)= yi – условие интерполяции Через заданные точки можно провести бесконечно много кривых, для каждой из которых выполнены все условия интерполяции. y1 x1 y2 x2 … … yi xi … … yn xn
Глобальная интерполяция Функция f(x) интерполируется на всем интервале [a; b] обычно с помощью единого интерполяционного полинома. Если количество узловых точек равно n, то степень полинома равна n-1, т.е. на единицу меньше количества узлов интерполяции.
Интерполяция полиномом Фактически задача сводится к определению коэффициентов интерполяционного полинома на основании значений функции в узловых точках. Существуют различные численные методы решения данной задачи. В пакете MATLAB имеются специальные функции для построения интерполяционного полинома.
Визуализация интерполируемой функции в MATLAB Для примера возьмем исходную функцию Построим график данной функции: >> x=linspace(2,7,1000); >> y=sin(2.*x).*sin(x); >> plot(x,y)
Интерполяция полиномом: постановка задачи Пусть заданы узлы интерполяции: xнач=2; xконечн=7; шаг Δx=0.5 Значения функции f(x) определяются в узлах интерполяции по формуле y=sin(2x)*sin(x). Требуется найти интерполяционный полином, график которого проходит через узловые точки интерполяции. Чтобы найти этот полином, надо определить его коэффициенты ai. Для определения коэффициентов, составляется система уравнений. Если количество узлов интерполяции 11 (в нашем случае), то полином 10 степени. Количество коэффициентов аi также равно11.
Интерполяция полиномом: расчет коэффициентов матричным методом Система уравнений для нахождения переменных ai в матричном виде: W*A=Y, где A- вектор-столбец переменных ai ( коэффициентов полинома): А=(a1 ; a2 ; a3 ; a4; a5 ; a6 ; a7; a8 ; a9 ; a10 ; a11). W=матрица Вандермонда Wij=xin-j степени аргументов в узлах интерполяции Y- вектор-столбец значений f(x) в узлах интерполяции Y=(y1 ;y2 ;y3 ;y4 ;y5 ;y6; y7 ;y8 ;y9 ;y10 ;y11) Матричное решение системы линейных алгебраических уравнений: A=inv(W)*Y W=
Блок-схема программы интерполяции полиномом
Расчет коэффициентов полинома матричным методом в MATLAB Решаем данную систему в MATLAB: >> x=2:0.5:7; % задаем узловые точки >> y=sin(2*x).*sin(x); % определяем значения y в узловых >> % точках >> W=vander(x); % определяем матрицу Вандермонда >> y1=y'; % преобразуем вектор-строку y в вектор->>%столбец >> A=inv(W)*y1; % решаем систему уравнений, находим >>% значения коэффициентов интерполяционного полинома
Значения коэффициентов полинома A = 0.00 -0.08 1.39 -14.24 91.86 -386.76 1067.07 -1882.63 1993.82 -1100.64 215.23 a1=0; a2=-0.08; a3=1.39; a4=-14.24; a5=91.86; a6=-386.76; a7=1067.07; a8=-1882.63; a9=1993.82; a10=-1100.64; a11=215.23 Интерполяционный полином: Р=0*а^10-0.08*а^9+1.39а8- -14.24а7+91.86а6- -386.76*а^5+1067.07*а^4- -1882.63*а^3+199.82*а^2- -1100.64*а+215.23
Расчет промежуточных значений функции в MATLAB при помощи полинома >> p=a; % задаем полином с коэффициентами ai >> x1=2:0.25:7; % задаем промежуточные значения х >> y1=polyval(p,x1); % вычисляем значения полинома в >>%промежуточных х >> plot(x,y,x1,y1); %распечатываем графики табличной >> %функции и интерполирующей функции
Глобальная или локальная? Пример: >> x=2:0.25:15; % задаем узловые точки >> length(x)% определяем размерность массива х ans = 53 % массив х состоит из 53-х элементов >> y=sin(2*x).*sin(x); % задаем значения функции в узловых точках >> W=vander(x); % определяем матрицу Вандермонда >> y1=y‘ % преобразуем вектор-строку y в вектор столбец >> a=inv(W)*y1; % вычисляем коэффициенты Warning: Matrix is close to singular or badly scaled.% предупреждение Results may be inaccurate. RCOND = 2.790744e-074. Система плохо обусловлена. Возможна ошибка. Вывод: При большом количестве узлов интерполяции трудно вычислить интерполяционный полином. MATLAB выдает сообщение о некорректной постановке задачи. В этом случае применяется локальная интерполяция.
Линейная интерполяция Узловые точки соединяются отрезками прямых xi xi+1 Интерполяция сплайнами (spline (англ.) – гибкая линейка) Интерполяция квадратичными сплайнами - через узловые точки проводят отрезки квадратичной параболы. Интерполяция кубическими сплайнами - через узловые точки проводят отрезки кубической параболы. xi xi+1 xi-1 Локальная (кусочно-полиномиальная) интерполяция На каждом интервале [xi , xi+1 ] строится отдельный интерполяционный полином невысокой степени
Локальная интерполяция функции в MATLAB Для одномерной табличной интерполяции используется функция interp1: yi = interp1(x, y, xi, method) — позволяет с помощью параметра method задать метод интерполяции: 'nearest' — ступенчатая интерполяция (по соседним элементам); 'linear' — линейная интерполяция (принята по умолчанию); 'spline' — кубическая сплайн-интерполяция; 'cubic' или 'pchip' — интерполяция многочленами Эрмита. Выходным аргументом interp1 является вектор промежуточных значений аргументов.
Блок-схема программы для кусочно-полиномиальной интерполяции
Ступенчатая интерполяция и ее визуализация в MATLAB Используем встроенную функцию: «ступенчатая интерполяция» - интерполяция по соседним элементам. >> x=2:0.5:7; >> y=sin(2.*x).*sin(x); >> X=linspace(2,7,21); >>Y1=interp1(x,y,X,'nearest'); >> plot(x,y) >> hold on >> plot(X,Y1) . Значение в каждой проме-жуточной точке принимается рав-ным ближайшему значению, заданному в таблице.
Линейная интерполяция и ее визуализация в MATLAB >> Y=interp1(x,y,X,'linear'); >> plot(x,y) >> hold on >> plot(X,Y) Используем встроенную функцию: «линейная интерполяция». Соединение соседних точек отрезками прямых - табличные данные приближаются ломаной линией .
Интерполяция кубическими сплайнами и ее визуализация в MATLAB >> Y2=interp1(x,y,X,'spline'); >> plot(x,y) >> hold on >> plot(X,Y2) Интерполяция кубическими сплайнами – кривыми 3ей степени - обеспечивает отсутствие сильных перегибов кривых в узлах интерполяции. Непрерывность функции и ее 1 и 2 производных на всем интервале интерполирования. Используем встроенную функцию: «интерполяция кубическими сплайнами».
Интерполяция многочленами Эрмитами ее визуализация в MATLAB >> Y3=interp1(x,y,X,'cubic'); >> plot(x,y) >> hold on >> plot(X,Y3) Используем встроенную функцию: «многочленами Эрмита». Аналогично интерполированию куби-ческими сплайнами , только в узлах интерполяции должны быть определены значения 1 и 2 производных интерполирую-щих функций.
Сравнение погрешностей методов интерполяции функции yi=sin(2*X).*sin(X) в MATLAB Файл-функция для определения погрешности: function p=pogr(X,Y); yi=sin(2.*X).*sin(X); P=abs(yi-Y); p=P; max (p) end
Погрешность ступенчатой интерполяции: >> X=linspace(2,7,51); >> Y=interp1(x,y,X,'nearest'); >> p=pogr(X,Y) ans = 0.3752
Погрешность линейной интерполяции: >> X=linspace(2,7,51); >> Y1=interp1(x,y,X,'linear'); >> p1=pogr(X,Y1) ans = 0.1400
Погрешность интерполяции кубическими сплайнами: >> X=linspace(2,7,51); >> Y2=interp1(x,y,X,'spline'); >> p2=pogr(X,Y2) ans = 0.0545
Погрешность интерполяции многочленами Эрмита: >> X=linspace(2,7,51); >> Y3=interp1(x,y,X,'cubic'); >> p3=pogr(X,Y3) ans = 0.1046 Вывод: наименьшую погрешность дает интерполяция кубическими сплайнами
Аппроксимация функций Аппроксимацией (приближением) функции f(x) называется нахождение такой функции g(x) (аппроксимирующей функции), которая была бы близка к заданной. Особенность метода - аппроксимирующая функция может быть произвольной. Наиболее часто встречаются аппроксимация прямой линией (линейная регрессия), аппроксимация полиномом (полиномиальная регрессия), аппроксимация линейной комбинацией произвольных функций.
Аппроксимация функций в MATLAB Одна из наиболее известных аппроксимаций — полиномиальная. В системе MATLAB определены функции аппроксимации данных полиномами по методу наименьших квадратов — полиномиальной регрессии. Это выполняет функция, приведенная ниже: – polyfit(x, y, n) – возвращает вектор коэффициентов полинома р(х) степени n, который с наименьшей среднеквадратичной погрешностью аппроксимирует функцию у(х).
Аппроксимация функции y=sin(2x)·sin(x) полиномами >> x=2:0.5:7; >> y=sin(2.*x).*sin(x); >> p=polyfit(x,y,3) % аппроксимация полиномом 3-ей степени p = -0.0113 0.1514 -0.3862 -0.3429 >> x1=2:0.1:7; >> y1=polyval(p,x1); >> plot(x,y,x1,y1)
Аппроксимация полиномом 3-ей степени
Аппроксимация полиномом 5-й степени >> p1=polyfit(x,y,5) % аппроксимация полиномом 5-й степени p1 = 0.0444 -1.0084 8.7876 -36.5622 72.4682 -55.0805 >> y2=polyval(p1,x1); >> plot(x,y,x1,y2)
Аппроксимация полиномом 5-й степени
Аппроксимация полиномом 7-й степени
Аппроксимация полиномом 10-й степени
Влияние степени полинома на вид приближающей функции При увеличении степени полинома график приближающей кривой лучше отражает форму графика аппроксимирующей функции. Если степень полинома равна n-1, где n – количество узловых точек, то аппроксимирующая кривая проходит через узловые точки, т.е. аппроксимация превращается в интерполяцию. Если степень полинома больше n, то задача имеет бесконечное множество решений, и MATLAB выводит сообщение об этом. При этом форма приближающей кривой на границах интервала далека от формы кривой аппроксимируемой функции.
Аппроксимация полиномом 20-й степени
Метод наименьших квадратов Для оценки «близости» функций выбирают тот или иной критерий согласия. Эти критерии основаны на использовании той или иной метрики, т.е. способа введения расстояния между функциями. Критерием близости в методе наименьших квадратов является требование минимальности суммы квадратов отклонений значений аппроксимирующей функции от экспериментальных точек:
Нахождение приближающей функции в виде линейной функции методом наименьших квадратов Пусть функция f(x) в точках xi принимает значения yi. Требуется найти аппроксимирующую (приближающую) функцию g(x). Приближающую функцию будем определять в виде: g(x)=a·x+b Для нахождения данной функции, необходимо определить параметры a и b.
Нахождение приближающей функции в виде линейной функции методом наименьших квадратов В результате математических выкладок по нахождению экстремума этой функции получаем систему уравнений , с помощью которой определяем параметры а и b линейной приближающей функции g(x)=a·x+b. Для применения критерия близости в методе наименьших квадратов необходимо найти минимум функции:
Система уравнений Мх, Мх2, Мхy – коэффициенты уравнений, определяемые по формулам:
Задача В одной системе координат построить графики табличной функции и аппроксимирующей функции. Решим задачу матричным методом в пакете MATLAB
Блок-схема программы нахождения приближающей функции в виде линейной функции методом наименьших квадратов
Программа для нахождения приближающей функции в виде линейной функции методом наименьших квадратов M-файл нахождения приближающей линейной функции: clear all % очищаем рабочую область x=input(‘вектор значений xi>>');% диалоговый ввод переменных хi y=input(‘вектор значений yi>>'); % диалоговый ввод переменных yi n=length(x) % определение размерности вектора х Mx=sum(x)/n % вычисление коэффициента системы уравнений My=sum(y)/n % вычисление коэффициента системы уравнений xy=x.*y; % поэлементное перемножение векторов x и y Mxy=sum(xy)/n % вычисление коэффициента системы уравнений x2=x.^2; % поэлементное возведение в квадрат вектора х Mx2=sum(x2)/n % вычисление коэффициента системы уравнений A=[Mx2 Mx;Mx 1]; % введение матрицы коэффициентов левой части системы… уравнений B=[Mxy;My]; % введение матрицы коэффициентов правой части системы уравнений ab=inv(A)*B; % решение системы уравнений a=ab(1) % угловой коэффициент приближающей линейной функции b=ab(2) % свободный член приближающей линейной функции Y=a*x+b; % аналитическое выражение для приближающей линейной функции plot(x,y,x,Y) % построение графиков исходной и приближающей функций
Нахождение приближающей функции в виде линейной функции (задача) вектор значений хi>>[1 1.5 2.5 3.0 4.5 5.1 6.2] вектор значений yi>>[67 101 168 202 301 334 404] n = 7 Mx = 3.400000000000000 My = 2.252857142857143e+002 Mxy = 9.724571428571428e+002 Mx2 = 14.742857142857142 a = 64.874326750448859 b = 4.713003334187988
Визуализация линейной приближающей функции
Путем построения кубического сплайна с помощью стандартной MATLAB процедуры interp1 определить значение термо-ЭДС термоэлектрического преобразователя ПП(S) при температуре 180оС, если известны значения термо-ЭДС при следующих температурах: Задача 1 на закрепление новой темы
Решение задачи 1 >> t=50:50:200; >> E=[0.299 0.645 1.029 1.44]; >> t1=180; >> E1=interp1(t,E,t1,'spline') E1 = 1.2731
С помощью стандартных MATLAB процедур polyfit и polyval определить сопротивление платинового термопреобразователя сопротивления при температуре 780оС путем аппроксимации полиномом 2-го порядка статической характеристики, которая представлена в таблице: Задача 2 на закрепление новой темы
>> t=650:50:950; >> R=[166.55 174.46 182.23 189.86 197.33 204.66 211.85]; >> p=polyfit(t,R,2) p = -0.0000 0.1976 50.4029 >> R1=polyval(p,780) R1 = 186.5429 Решение задачи 2