Этот пример показывает, что экстраполирование данных с помощью полиномов даже небольшой степени опасно и ненадежно.
Этот пример является более старым, чем MATLAB®. Это запустилось как упражнение в Компьютерных Методах для Математических вычислений, Форсайтом, Малкольмом и Молером, опубликованным Prentice Hall в 1 977.
Теперь, MATLAB делает намного легче отличаться параметры и видеть результаты, но базовые математические принципы неизменны.
Создайте и постройте два вектора с данными о переписи США от 1 910 до 2000.
% Time interval t = (1910:10:2000)'; % Population p = [91.972 105.711 123.203 131.669 150.697... 179.323 203.212 226.505 249.633 281.422]'; % Plot plot(t,p,'bo'); axis([1910 2020 0 400]); title('Population of the U.S. 1910-2000'); ylabel('Millions');
Каково ваше предположение для генеральной совокупности в году 2010?
p
p = 10×1
91.9720
105.7110
123.2030
131.6690
150.6970
179.3230
203.2120
226.5050
249.6330
281.4220
Соответствуйте данным полиномом в t
и используйте его, чтобы экстраполировать генеральную совокупность в t = 2010
. Получите коэффициенты в полиноме путем решения линейной системы уравнений, включающих 11 11 матрицу Вандермонда, с элементами как степени масштабированного времени, A(i,j) = s(i)^(n-j)
.
n = length(t); s = (t-1950)/50; A = zeros(n); A(:,end) = 1; for j = n-1:-1:1 A(:,j) = s .* A(:,j+1); end
Получите коэффициенты c
для полинома степени d
, который соответствует данным p
путем решения линейной системы уравнений, включающих последние столбцы d+1
матрицы Вандермонда:
A(:,n-d:n)*c ~= p
Если d < 10
, то больше уравнений, чем неизвестные существует, и решение методом наименьших квадратов, является соответствующим.
Если d == 10
, то можно решить уравнения точно и полином на самом деле, интерполирует данные.
В любом случае используйте оператор наклонной черты влево, чтобы решить систему. Коэффициенты для кубического соответствия:
c = A(:,n-3:n)\p
c = 4×1
-5.7042
27.9064
103.1528
155.1017
Теперь оцените полином в каждый год от 1 910 до 2010 и постройте результаты.
v = (1910:2020)'; x = (v-1950)/50; w = (2010-1950)/50; y = polyval(c,x); z = polyval(c,w); hold on plot(v,y,'k-'); plot(2010,z,'ks'); text(2010,z+15,num2str(z)); hold off
Сравните кубическое соответствие с биквадратным Уведомлением, что экстраполируемая точка очень отличается.
c = A(:,n-4:n)\p; y = polyval(c,x); z = polyval(c,w); hold on plot(v,y,'k-'); plot(2010,z,'ks'); text(2010,z-15,num2str(z)); hold off
Когда степень увеличивается, экстраполяция становится еще более ошибочной.
cla plot(t,p,'bo') hold on axis([1910 2020 0 400]) colors = hsv(8); labels = {'data'}; for d = 1:8 [Q,R] = qr(A(:,n-d:n)); R = R(1:d+1,:); Q = Q(:,1:d+1); c = R\(Q'*p); % Same as c = A(:,n-d:n)\p; y = polyval(c,x); z = polyval(c,11); plot(v,y,'color',colors(d,:)); labels{end+1} = ['degree = ' int2str(d)]; end legend(labels, 'Location', 'NorthWest') hold off