Подгонка кривой 5 баллов

Я пытаюсь подогнать кривую под 5 точек в C. Я использовал этот код из предыдущего сообщения (Может ли кто-нибудь упростить мне это уравнение?) сделать 4 балла, но теперь мне нужно добавить еще один балл.

// Input data: arrays x[] and y[]
// x[1],x[2],x[3],x[4] - X values
// y[1],y[2],y[3],y[4] - Y values

// Calculations
A = 0
B = 0
C = 0
D = 0
S1 = x[1] + x[2] + x[3] + x[4]
S2 = x[1]*x[2] + x[1]*x[3] + x[1]*x[4] + x[2]*x[3] + x[2]*x[4] + x[3]*x[4]
S3 = x[1]*x[2]*x[3] + x[1]*x[2]*x[4] + x[1]*x[3]*x[4] + x[2]*x[3]*x[4]
for i = 1 to 4 loop
   C0 = y[i]/(((4*x[i]-3*S1)*x[i]+2*S2)*x[i]-S3)
   C1 = C0*(S1 - x[i])
   C2 = S2*C0 - C1*x[i]
   C3 = S3*C0 - C2*x[i]
   A = A + C0
   B = B - C1
   C = C + C2
   D = D - C3
end-loop

// Result: A, B, C, D

Я пытался скрыть это до 5-точечной кривой, но у меня возникли проблемы с выяснением того, что происходит внутри цикла:

// Input data: arrays x[] and y[]
// x[1],x[2],x[3],x[4],x[5] - X values
// y[1],y[2],y[3],y[4],y[5] - Y values

// Calculations
A = 0
B = 0
C = 0
D = 0
E = 0
S1 = x[1] + x[2] + x[3] + x[4]
S2 = x[1]*x[2] + x[1]*x[3] + x[1]*x[4] + x[2]*x[3] + x[2]*x[4] + x[3]*x[4]
S3 = x[1]*x[2]*x[3] + x[1]*x[2]*x[4] + x[1]*x[3]*x[4] + x[2]*x[3]*x[4]
S4 = x[1]*x[2]*x[3]*x[4] + x[1]*x[2]*x[3]*[5] + x[1]*x[2]*x[4]*[5] + x[1]*x[3]*x[4]*[5] + x[2]*x[3]*x[4]*[5]

for i = 1 to 4 loop
   C0 = ??
   C1 = ??
   C2 = ??
   C3 = ??
   C4 = ??
   A = A + C0
   B = B - C1
   C = C + C2
   D = D - C3
   E = E + C4
end-loop

// Result: A, B, C, D, E

любая помощь в заполнении C0...C4 будет оценена по достоинству. Я знаю, что это связано с матрицами, но я не смог понять это. примеры с псевдокодом или реальным кодом были бы наиболее полезными.

Спасибо


person user3550036    schedule 19.05.2014    source источник
comment
Равноудалены ли значения x?   -  person John Alexiou    schedule 19.05.2014
comment
Спасибо за ответ. да, они более или менее расположены на одинаковом расстоянии друг от друга, окончательная кривая должна выглядеть примерно как колоколообразная кривая с некоторыми вариациями, конечно.   -  person user3550036    schedule 19.05.2014
comment
@ user3550036 К чему именно вы приближаетесь? Вы говорите, что окончательная кривая должна выглядеть как кривая колокола — если вы имеете в виду, что рассматриваемая функция является гауссовой (то есть нормальной) плотностью, это имеет большое значение. Если это так, перейдите непосредственно к гауссовой плотности вместо полинома — полиномы не обязательно будут положительными везде, даже если все точки данных положительны.   -  person Robert Dodier    schedule 19.05.2014


Ответы (2)


Я не хочу упускать эту возможность обобщить. :)

Вместо этого мы немного познакомимся с полиномами Лагранжа и Метод разделенной разности Ньютона их вычисления.

Многочлены Лагранжа

Учитывая n+1 точек данных, интерполирующий полином равен

Интерполяционное полиномиальное определение

где l_j(i)

l_i определение.

Это означает, что мы можем найти полином, приближающий n+1 точек, независимо от интервала и т. д., просто суммируя эти полиномы. Однако это немного сложно, и я бы не хотел делать это на C. Давайте посмотрим на Многочлены Ньютона.

Многочлены Ньютона

То же самое начало, учитывая n+1 точки данных, аппроксимирующий полином будет

введите здесь описание изображениягде каждый n(x)

введите здесь описание изображенияс коэффициентом

введите здесь описание изображения, являющееся разделенной разницей.

Окончательная форма выглядит так

введите здесь описание изображения.

Как видите, формула довольно проста, учитывая значения разделенной разницы. Вы просто делаете каждую новую разделенную разницу и умножаете на каждую точку до сих пор. Следует отметить, что вы получите полином степени n из n+1 точек.

Разделенная разница

Все, что осталось, это определить разделенную разницу, которую действительно лучше всего объясняют эти две картинки:

введите здесь описание изображения

а также

введите здесь описание изображения.

С этой информацией реализация C должна быть разумной. Я надеюсь, что это поможет, и я надеюсь, что вы узнали что-то! :)

person Logan    schedule 19.05.2014
comment
Привет, спасибо за вашу помощь. В ваших элементах, обведенных красным, я не вижу y0, но в приведенной выше окончательной форме уравнения используется y0. не могли бы вы уточнить, как эти [y0], [y0,y1] переводятся в то, что у вас выделено красным? - person user3550036; 19.05.2014
comment
Это просто разница в индексации. Один с красными кружками начинает индексацию с 1, а другие начинают с 0. :) - person Logan; 19.05.2014
comment
[y0] было бы y1, а [y0,y1] было бы y2,1, если бы вы перевели нотацию. - person Logan; 19.05.2014
comment
А если вы хотите получить колоколообразную кривую, используйте полиномиальную интерполяцию для данных x[k], sqrt(C-log(y[k])), чтобы получить полином p(x). Тогда интерполирующая кривая будет exp(C-p(x)^2). Выберите C больше, чем, например, вдвое или плюс один максимальное значение любого log(y[k]). - person Lutz Lehmann; 19.05.2014

Если значения x расположены на одинаковом расстоянии друг от друга с x2-x1=h, x3-x2=h, x4-x3=h и x5-x4=h, то

C0 = y1;
C1 = -(25*y1-48*y2+36*y3-16*y4+3*y5)/(12*h);
C2 = (35*y1-104*y2+114*y3-56*y4+11*y5)/(24*h*h);
C3 = -(5*y1-18*y2+24*y3-14*y4+3*y5)/(12*h*h*h);
C4 = (y1-4*y2+6*y3-4*y4+y5)/(24*h*h*h*h);

y(x) = C0+C1*(x-x1)+C2*(x-x1)^2+C3*(x-x1)^3+C4*(x-x1)^4
// where `^` denotes exponentiation (and not XOR).
person John Alexiou    schedule 19.05.2014
comment
Обратите внимание, что приведенный выше полином 4-го порядка, и вместо этого вам могут понадобиться 4 кубических сплайна 3-го порядка. Исходный вопрос не уточняется. - person John Alexiou; 19.05.2014
comment
Спасибо за помощь. Мне нужен полином 4-го порядка. Что бы я сделал, чтобы заменить hs, если они не расположены точно на одинаковом расстоянии? не могли бы вы просто заменить их на h1, h2, h3, h4, соответствующие (x2-x1) = h1, (x3-x2) = h2, x4-x3 = h3, x5-x4 = h4? или это математически не разрешено - person user3550036; 19.05.2014
comment
Если они не разделены точно на h, то инверсия матрицы 4-го порядка матрицы Вандермонда далека. более сложный. - person John Alexiou; 19.05.2014