Как воссоздать функцию построенную сплайном матлаб
Для реализации сплайнов в Matlab есть встроенная функция interpl следующего вида:
yi = interpl ( x , y , xi , method ), где x – массив абсцисс экспериментальных точек, y – массив ординат экспериментальных точек, xi – точки, в которых необходимо вычислить значение с помощью сплайна, method – определяет метод построения сплайна. Параметр method может принимать одно из следующих значений:
‘ linear ’ – линейная интерполяция;
‘ spline ’ – кубический сплайн;
‘ nearest ’ – интерполяция по соседним точкам – это метод построения кусочной функции, при котором значение в любой точке равно значению в ближайшей узловой точке.
Функция interpl ( x , y , xi , method ) возвращает значения интерполирующей функции в точках xi .
В программе представлено решение задачи 1 с помощью функции interpl :
x=[0.43 0.48 0.55 0.62 0.7 0.75];
y=[1.635997 1.73234 1.87686 2.03345 2.35973];
t=[0.702 0.512 0.608];
ytn= interpl(x,y,t, ‘nearest’);
ytl= interpl(x,y,t, ‘linear’);
yts= interpl(x,y,t, ‘spline’);
Пример сплайн – интерполяции.
Пусть функция y = sin ( x ) задано таблично в 8 точках на отрезке [0,2 pi ]. Выполнить сплайн-интерполяцию данной функции.
% 1. Задание табличных значений интерполируемой функции в точках x ( i )
>> dx=2*pi/(N-1); % шаг дискретизации в таблице
% 2. Задание значения абсцисс точек, в которых вычисляется значение интерполяционного полинома
% 3. Вычисление интерполируемых значений функции Y в узлах сетки X
% 4. Визуализация результатов интерполяции
>> plot(x,y, ‘o’, X, Y); xlabel (‘X’); ylabel(‘Y’); grid on;
Вычисление погрешности между точным и интерполируемыми значениями
>> Y _ t = sin ( X ); % вычисление точных значений интерполируемой функции
% построение нового графика с белым фоном
>> figurel = figure (‘Color’, [1 1 1]);
>> figurel = figure (‘Color’, [1 1 1]);
>> plot (X, abc (Y_t-Y)); xlabel (‘X’); ylabel (‘Y’);
>> title (‘Абсолютная погрешность’);
В результате получим график функции
Погрешность на краях отрезка существенно выше, чем в середине отрезка. Это происходит потому, что в функции spline , используемой в данном примере, не учитываются граничные условия .
Как воссоздать функцию построенную сплайном матлаб
>> x = [1 2 3 4 5 6 7 8];
>> y = [-1.1 0.2 0.5 0.8 0.7 0.6 0.4 0.1];
ans = 0.1143 -0.2393
y ( x ) = 0.1143 x – 0.2393
>> plot ( x , y ), grid
>> x =-3* pi : pi /100:3* pi ;
>> x = [1 2 3 4 5 6 7 8];
>> y = [-1.1 0.2 0.5 0.8 0.7 0.6 0.4 0.1];
>> x= [1 3 2 9 6 8 4 9 -6];
function y = f1(x, a)
function y = parab (x)
y = rectpuls ( t , w ).
>> y = 0.75*rectpuls(t-3,2) + 0.5*rectpuls(t-8, 0.4)…
y= tripuls (t, w, s)
>> y = 0.75*tripuls(t-1, 0.5) + 0.5*tripuls(t-5, 0.5, -1)…
+ 1.35*tripuls(t-3, 0.8, 1);
y= gauspuls (t, fc, bw),
>> y = 0.75*gauspuls(t-3, 1, 0.5);
y = square ( t , duty ),
y = sawtooth ( t , width )
>> y = 0.7* sawtooth (pi*t/5, 0.5);
h = freqs ( b , a , w )
>>a1(1) = 1; a1(2) = 2*dz*om0; a1(3) = om0^2; b1(1)=A;
y = fft ( x , n ); x = ifft ( y , n ),
>> f 1 = -500 : 0.5 :500;
>> plot (f1(970:1030), a(970:1030)); grid
>> plot (f1(970:1030), a(970:1030)); grid
y = filter (b, a, x),
y(k) = b(1)*x(k) + b(2)*x(k-1) + b(nb+1)*x(k-nb)-
>> T2 = 0.2; A2 = 10; eps = pi/4;
>> T1 = 1; Tf = T1; dz = 0.05;
>> om0=2*pi/Tf; A=1; oms = om0*Ts;
>> plot(t(1002:end), y(t(1002:end),t(1002:end), Yp(1002:end), grid,…
>> [C,f] = psd (y1, dovg, Fmax);
>> psd ( y 1, dovg , Fmax )
>> tau = -10+Ts : Ts : 10;
>> s1r = round (length ( R) /2 – lt/2;
>> s2r = round (length ( R) /2 + lt/2 -1;
>> plot (tau, R (s1r:s2r), grid
>> v = tf ([1 4], [1 2 100])
1. Гультяев А, Имитационное моделирование в среде Windows. – С-Пб.: “КОРОНАпринт”, 1999. 2. Лазарев Ю. MatLAB 5.x. – Киев: «Ирина», BHV, 2000.
6. Discrete Filter.
7. Discrete Transfer Fen.
8. Discrete Zero-Pole.
(Initial condition source).
4. Zero — Pole .
§ I = trapz (x, y)
x = [-0.5 -0.4677 -0.4343 -0.3952 -0.3493 -0.2951 -0.2278 -0.1431 -0.0394 0.0803 0.2099 0.3445 0.4816 0.6225 0.7647 0.9009 1];
y = [16 14.6264 13.2183 11.8263 10.4538 9.1135 7.8345 6.6780 5.7083 4.9382 4.3355 3.8454 3.4550 3.1773 2.9703 2.7792 2.7];
Сплайны
Английское слово spline можно перевести как «гибкая линейка». Сплайны, в отличие полиномиальной интерполяции предыдущих пунктов, представляют собой полиномы невысокой степени, например первой или третьей. Строятся они нс глобально для всего отрезка интерполяции, а применительно к более мелкому отрезку между парой узлов интерполяции. При этом принято обеспечивать совпадение значений соседних сплайнов в узлах интерполяции и, быть может, непрерывность первых или вторых производных сплайнов интерполирующей функции. Техника сплайнов может быть еще охарактеризована как процедура построения кусочно-полиномиальной интерполяции.
Более подробно рассмотрим наиболее употребляемый кубический сплайн. Пусть искомая функция / задана в п узлах интерполяции Хц . Х„ своими значениями yi, . у„. На отрезке [x,,Xj+i] определим кубический сплайн вида Sj(x) = а, + Ь,Х + CjX 2 + d[X 3 , Всего, таким образом, можно определить П-1
сплайнов S,(х), где ( = 1, . п-1. Для определения 4(n—1) неизвестных а,, Ь„ С, d„ i = = 1, . П-1, составим П условий совпадения значений сплайна со значениями интерполируемой функции:
П-2 условий непрерывности сплайна: 
П-2 условий непрерывности производной сплайна:
П-2 условий непрерывности второй производной сплайна:
Всего, таким образом, на 4л-4 неизвестных сплайна накладывается 4Л-6 уравнений. Недостающие два уравнения обычно определяют, исходя из тех или иных граничных условий. Предположим, например, что функция / на своих границах удовлетворяет условиям /»(х„) = 0, тогда имеем
следующую недостающую пару уравнений:
Уравнения (2.9) — (2.13) однозначно определяют кубический сплайн S=s(x). Для изучения процесса сходимости сплайна к интерполируемой функции сформулируем без доказательства теорему. Для формулировки теоремы введем ряд предположений.
Рис. 2.4. Пример построения кубического сплайна по 41 узлу интерполяции
Пусть интерполируемая функция задана на равномерной сетке Х„ 1 = 1, . П, X, = (l-l)/l, где h = l/(n—1) — шаг сетки. Пусть функция f имеет четвертую производную, т. е. относится к классу функций С (4) и, кроме того,
Г(0) = /»(1) = О . Пусть кубический сплайн s(x) удовлетворяет тем же условиями, т. е. s”(0) = S»(l) = 0, тогда после введения обозначений
можно сформулировать теорему об оценке ошибки интерполяции для функции Дх) и ее производных f (х), Д'(х).
Теорема [3, с. 144]. Для f € С (4) справедливы оценки
Из оценок (2.14) следует, что при h —> 0 (п —> со) нс только сплайн сводится к интерполируемой функции, но и его первая и вторая производный сходятся к соответствующим производным интерполируемой функции.
Займемся теперь программированием сплайнов. В листинге 2.4 приведен код программы для построения кубического сплайна. На рис. 2.4 приведен пример расчета сплайна.
%Интерполяция с помощью сплайна %определим вектор узлов интерполяции х =0: 0. 0 2 5: 1;
%определим значения интерполируемой функции,
%считая эти значения случайными величинами, %распределенными по нормальному закону У=[ 1 ;
for i =1:Iengt h(х) у =[ у г a n d п ]; end
%вводим сетку на отрезке интерполяции
%обращаемся к стандартной процедуре MATLAB yv=i nterpl(x,y,xv,’spl i ne’ );
piot(x,у, 1 *’,xv,yv);
Стандартная функция interpl в MATLAB дает возможность также построить кусочно-постоянные сплайны полиномами нулевой степени (способ ‘nearest’) и линейные сплайны полиномами первой степени (способ linear 1 ). На рис. 2.5, а приведен пример кусочно-постоянного сплайна, а на рис. 2.5, б пример линейного сплайна.
Рис. 2.5. Сплайн:
а — кусочно-постоянный; б — линейный
Функция interpl из листинга 4 допускает еще два способа построения сплайнов: ‘pchip’ и ‘cubic’. В определенном смысле эти способы дают более высокого качества сплайны, чем предыдущие процедуры. Предлагается студентам разобраться с этими режимами построения сплайнов самостоятельно.
Как воссоздать функцию построенную сплайном матлаб
Suppose you want to interpolate some smooth data, e.g., to
You can use the cubic spline interpolant obtained by
and plot the spline, along with the data, with the following code:
This produces a figure like the following.
Cubic Spline Interpolant of Smooth Data
This is, more precisely, the cubic spline interpolant with the not-a-knot end conditions, meaning that it is the unique piecewise cubic polynomial with two continuous derivatives with breaks at all interior data sites except for the leftmost and the rightmost one. It is the same interpolant as produced by the MATLAB ® spline command, spline(x,y) .
Periodic Data
The sine function is 2π-periodic. To check how well your interpolant does on that score, compute, e.g., the difference in the value of its first derivative at the two endpoints,
which is not so good. If you prefer to get an interpolant whose first and second derivatives at the two endpoints, 0 and 4*pi , match, use instead the command csape which permits specification of many different kinds of end conditions, including periodic end conditions. So, use instead
for which you get
Output is ans = 0 as the difference of end slopes. Even the difference in end second derivatives is small:
Output is ans = -4.6074e-015 .
Other End Conditions
Other end conditions can be handled as well. For example,
provides the cubic spline interpolant with breaks at the and with its slope at the leftmost data site equal to 3, and its second derivative at the rightmost data site equal to -4.
General Spline Interpolation
If you want to interpolate at sites other than the breaks and/or by splines other than cubic splines with simple knots, then you use the spapi command. In its simplest form, you would say sp = spapi(k,x,y) ; in which the first argument, k , specifies the order of the interpolating spline; this is the number of coefficients in each polynomial piece, i.e., 1 more than the nominal degree of its polynomial pieces. For example, the next figure shows a linear, a quadratic, and a quartic spline interpolant to your data, as obtained by the statements
Spline Interpolants of Various Orders of Smooth Data
Even the cubic spline interpolant obtained from spapi is different from the one provided by csapi and spline . To emphasize their difference, compute and plot their second derivatives, as follows:
This gives the following graph:
Second Derivative of Two Cubic Spline Interpolants of the Same Smooth Data
Since the second derivative of a cubic spline is a broken line, with vertices at the breaks of the spline, you can see clearly that csapi places breaks at the data sites, while spapi does not.
Knot Choices
It is, in fact, possible to specify explicitly just where the spline interpolant should have its breaks, using the command sp = spapi(knots,x,y) ; in which the sequence knots supplies, in a certain way, the breaks to be used. For example, recalling that you had chosen y to be sin(x) , the command
provides a cubic Hermite interpolant to the sine function, namely the piecewise cubic function, with breaks at all the x(i) ‘s, that matches the sine function in value and slope at all the x(i) ‘s. This makes the interpolant continuous with continuous first derivative but, in general, it has jumps across the breaks in its second derivative. Just how does this command know which part of the data value array [y cos(x)] supplies the values and which the slopes? Notice that the data site array here is given as [x x] , i.e., each data site appears twice. Also notice that y(i) is associated with the first occurrence of x(i) , and cos(x(i)) is associated with the second occurrence of x(i) . The data value associated with the first appearance of a data site is taken to be a function value; the data value associated with the second appearance is taken to be a slope. If there were a third appearance of that data site, the corresponding data value would be taken as the second derivative value to be matched at that site. See Constructing and Working with B-form Splines for a discussion of the command augknt used here to generate the appropriate «knot sequence».
Smoothing
What if the data are noisy? For example, suppose that the given values are
Then you might prefer to approximate instead. For example, you might try the cubic smoothing spline, obtained by the command