MATLAB - Попытка сделать полиномиальную интерполяцию, но я получаю матрицу, полную нулей (почти)
Я пытаюсь написать код, который выполняет полиномиальную интерполяцию, но я не могу заставить его работать. Я студент и следую логике из этого видео https://www.youtube.com/watch?v=VpI-wC94RKw, чтобы воссоздать его в виде кода в Matlab, но всякий раз, когда я пытаюсь создать свой Собственная версия матрицы, показанной в видео, вместо этого я получаю матрицу, почти исключительно заполненную нулями (кроме одного элемента). Я не могу понять, почему это происходит.
Мой код:
x=[150, 200, 300, 500, 1000, 2000, 99999]';
y=[2, 3, 4, 5, 6, 7, 8]';
function f = interPoly(x,y)
% Skapar en matris A var varje rad är [x_1^6, x_1^5,..., 1] till [x_n^6, x_n^5,..., 1]
A = [x.^6 x.^5 x.^4 x.^3 x.^2 x ones(numel(x),1) y];
% Gaussar matris A
R = rref(A);
% Plockar sista kolumnen ur R
c = R(:,end);
f = c(1)*x.^6+c(2)*x.^5+c(3)*x.^4+c(4)*x.^3+c(5)*x.^2+c(6)*x+c(7);
end
(Матрица "А" здесь проблематична. Функция, которую я получаю в конце, также просто заполнена нулями в качестве значений. Также извините за комментарии на шведском)
У меня есть 7 значений в x и y и, следовательно, полином 6-го порядка, но я не знаю, какой должна быть константа во втором последнем столбце, поэтому я просто поставил несколько из них (я новичок в этом, поэтому я немного неуверен в логике).
В любом случае я пробовал использовать ту же функцию с некоторыми другими входными данными, и она работала нормально.
Альтернативные входные данные:
x=[0, 0.5, 1, 1.5, 2, 2.99, 3]';
y=[0, 0.52, 1.09, 1.75, 2.45, 3.5, 4]';
Это дает мне нули, потому что элементы переполняются (например, 99999^6 - очень большое число)? Я не очень понимаю, что здесь происходит и почему он отлично работает с другим набором входных данных. Помогите?
Спасибо!
РЕДАКТИРОВАТЬ: весь смысл этой задачи (заданной моей школой) состоит в том, чтобы сравнить метод "наименьших квадратов" (для которого я также написал код, но не опубликовал) с методом полиномиальной интерполяции (тот, что в коде выше). Последнее значение в 'x' выше должно быть бесконечностью (f(inf)=8), поэтому я просто заменил его действительно большим числом, поэтому оно не распределяется "равномерно". Есть лучший способ сделать это?
1 ответ
В вашем примере вы получаете следующую матрицу для R:
1.00000 0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 -0.00000 1.6925e-05
0.00000 1.00000 0.00051 0.00000 0.00000 0.00000 0.00000 0.00000 7.1286e-05
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 5.4078e-04
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 6.9406e-03
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 2.2098e-01
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 7.0000e+00
0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 8.0000e+00
Если система уравнений, которую вы определяете, имеет одно уникальное решение, то у вас должна быть единичная матрица (с этим дополнительным столбцом справа). В вашем случае у вас есть (из первых двух строк) уравнения
0 == 5.4078e-04
0 == 6.9406e-03
0 == 2.2098e-01
0 == 7.0000e+00
0 == 8.0000e+00
которые указывают, что ваша система не допускает решения. Это связано с вашим выбором входных данных, которые теоретически могут иметь решение, но на практике числовые свойства очень плохие. И действительно, если мы вычислим число условий cond(A) мы получаем значение 4.3691e+06 что очень плохо для такой маленькой матрицы. Скорее всего это связано с неравномерным распространением вашего x ценности. (Вы пытаетесь решить систему Вандермонда.) Если вы распределяете их более "равномерно", все выглядит намного лучше.
Кроме того, ваш код выглядит хорошо, но вы, вероятно, хотите оценить функцию интерполяции в значениях, отличных от тех, которые используются для их определения. Кроме того, этот "точный" подход, как известно, не очень идеален в численном отношении, поэтому, возможно, стоит использовать более устойчивый алгоритм для решения системы, чем rrefнапример qr разложение. Я предлагаю следующие изменения:
function f = interPoly(x,y,xnew)
A = [x.^6 x.^5 x.^4 x.^3 x.^2 x ones(numel(x),1)];
[q,r] = qr(A);
c = r\(q' * y);
disp(cond(A));
f = c(1)*xnew.^6+c(2)*xnew.^5+c(3)*xnew.^4+c(4)*xnew.^3+c(5)*xnew.^2+c(6)*xnew+c(7);
end
x=[150, 200, 300, 500, 1000, 2000, 99999]';
y=[2, 3, 4, 5, 6, 7, 8]';
disp(interPoly(x,y,x+eps));