MATLAB - Попытка сделать полиномиальную интерполяцию, но я получаю матрицу, полную нулей (почти)

Question

Я пытаюсь написать код, который выполняет полиномиальную интерполяцию, но не могу заставить его работать.Я студент, и я следую логике из этого видео 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), поэтому я просто заменил его действительно большим числом, поэтому оно не распределяется «равномерно».Есть ли лучший способ сделать это?

Answer 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.(Вы пытаетесь решить систему Vandermonde .) Если вы распределяете их более "равномерно", все выглядит намного лучше.

Кроме того, ваш код выглядит нормально, но вы, вероятно, хотите оценить интерполяционную функцию в значениях, отличных от тех, которые используются для их определения.Кроме того, известно, что этот «точный» подход не очень идеален в численном отношении, поэтому, возможно, стоило бы использовать более устойчивый алгоритм для решения системы, чем 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));

Попробуйте онлайн!