Interpolácia pozostáva z odhadu hodnot (bodov) funkcie medzi zadanými hodnotami (x0, y0), (x1, y1), ... (xi-1, yi-1), (xi, yi), (xi+1, yi+1), ... (xn, yn). Hodnoty sa interpolujú vždy medzi dvojicou zadaných hôdnot, čize pre nejaký interval. Pričom do výpočtu hodnôt z daného intervalu vtupujú hodboty z predchádzajúceho aj nasledujúcého intervalu, tým sa dosiahne hladši a presnejši priebeh interpolovanej funkcie. Podrobné vysvetlenie a odvodzovanie výrazov pre spline interpoláciu je voľne dostupné na webe. Takže len pre informáciu sú nižsie uvedené už len odvodené vzorce, ktoré sa využivajú v programe. Cieľom je vypočitať koeficienty "a, b, c, d" pre jednotlivé intervaly a pomocou polynómu tretieho rádu získať interpolované body funkcie a teda priebeh interpolovanej funkcie.
Polynóm tretieho rádu pre výpočet bodov funkcie:
Na výpočet jednotlivých koeficientov je potrebné poznat prvú a druhú deriváciu:
Po úprave (úprava je nutná aby v zdrojovom kóde nenastalo delenie nulov, v prípade ak je y súradnica funkcie rovnaká pre rôzne x):
Pričom pre krajné (hraničné) body platí:
Druhé derivácie:
Koeficienty pre jednotlive intevaly zadaných hodnôt sa vypočítajú podľa vzorcov:
Ako príklad je uvedená interpolácia pre 6 intervalov, tzn. pre 7 zadaných súradníc (x,y):
Graf interpolácie - zadané body (červená farba), vypočítané body (modrá farba):
Zdrojový kód
public bool splineInterpolation(out string errorMessage, double[] xInArray, double[] yInArray, out double[] yOutArray, out double[,] koeffArray)
{
errorMessage = "";
if (xInArray.Length < 3)
{
errorMessage = "Array Lenght Error: Vstupne pole (polia) ma nespravnu dlzku.";
yOutArray = new double[0];
koeffArray = new double[0, 0];
return false;
}
int pointCount = Convert.ToInt32(xInArray[xInArray.Length - 1]) - Convert.ToInt32(xInArray[0]) + 1;
yOutArray = new double[pointCount];
double f1x0_deriv = 0;
double f1x1_deriv = 0;
double f2x0_deriv = 0;
double f2x1_deriv = 0;
int intervalLenght = xInArray.Length - 1;
koeffArray = new double[intervalLenght, 4];
try
{
if (xInArray.Length != yInArray.Length)
{
errorMessage = "Array Lenght Error: Vstupne pole (polia) ma nespravnu dlzku.";
return false;
}
//Vypocet koeficientov pre jednotlive intervaly
for (int iInterval = 0; iInterval < intervalLenght; iInterval++)
{
if (iInterval == 0)
{
f1x1_deriv = (2 * (yInArray[iInterval + 2] - yInArray[iInterval + 1])
* (yInArray[iInterval + 1] - yInArray[iInterval]))
/ ((xInArray[iInterval + 2] - xInArray[iInterval + 1])
* (yInArray[iInterval + 1] - yInArray[iInterval])
+ (xInArray[iInterval + 1] - xInArray[iInterval])
* (yInArray[iInterval + 2] - yInArray[iInterval + 1]));
f1x0_deriv = ((3 * (yInArray[iInterval + 1] - yInArray[iInterval]))
/ (2 * (xInArray[iInterval + 1] - xInArray[iInterval])))
- (f1x1_deriv / 2);
f2x0_deriv = ((-2 * (f1x1_deriv + 2 * f1x0_deriv))
/ (xInArray[iInterval + 1] - xInArray[iInterval]))
+ ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
/ Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
f2x1_deriv = ((2 * (2 * f1x1_deriv + f1x0_deriv))
/ (xInArray[iInterval + 1] - xInArray[iInterval]))
- ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
/ Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
}
else
{
if (iInterval == intervalLenght - 1)
{
f1x0_deriv = (2 * (yInArray[iInterval + 1] - yInArray[iInterval])
* (yInArray[iInterval] - yInArray[iInterval - 1]))
/ ((xInArray[iInterval + 1] - xInArray[iInterval])
* (yInArray[iInterval] - yInArray[iInterval - 1])
+ (xInArray[iInterval] - xInArray[iInterval - 1])
* (yInArray[iInterval + 1] - yInArray[iInterval]));
f1x1_deriv = ((3 * (yInArray[iInterval + 1] - yInArray[iInterval]))
/ (2 * (xInArray[iInterval + 1] - xInArray[iInterval])))
- (f1x0_deriv / 2);
f2x0_deriv = ((-2 * (f1x1_deriv + 2 * f1x0_deriv))
/ (xInArray[iInterval + 1] - xInArray[iInterval]))
+ ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
/ Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
f2x1_deriv = ((2 * (2 * f1x1_deriv + f1x0_deriv))
/ (xInArray[iInterval + 1] - xInArray[iInterval]))
- ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
/ Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
}
else
{
f1x1_deriv = (2 * (yInArray[iInterval + 2] - yInArray[iInterval + 1])
* (yInArray[iInterval + 1] - yInArray[iInterval]))
/ ((xInArray[iInterval + 2] - xInArray[iInterval + 1])
* (yInArray[iInterval + 1] - yInArray[iInterval])
+ (xInArray[iInterval + 1] - xInArray[iInterval])
* (yInArray[iInterval + 2] - yInArray[iInterval + 1]));
f1x0_deriv = (2 * (yInArray[iInterval + 1] - yInArray[iInterval])
* (yInArray[iInterval] - yInArray[iInterval - 1]))
/ ((xInArray[iInterval + 1] - xInArray[iInterval])
* (yInArray[iInterval] - yInArray[iInterval - 1])
+ (xInArray[iInterval] - xInArray[iInterval - 1])
* (yInArray[iInterval + 1] - yInArray[iInterval]));
f2x0_deriv = ((-2 * (f1x1_deriv + 2 * f1x0_deriv))
/ (xInArray[iInterval + 1] - xInArray[iInterval]))
+ ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
/ Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
f2x1_deriv = ((2 * (2 * f1x1_deriv + f1x0_deriv))
/ (xInArray[iInterval + 1] - xInArray[iInterval]))
- ((6 * (yInArray[iInterval + 1] - yInArray[iInterval]))
/ Math.Pow(xInArray[iInterval + 1] - xInArray[iInterval], 2));
}
}
//koeficienty a, b, c, d - stlpce, intervaly - riadky
koeffArray[iInterval, 0] = (f2x1_deriv - f2x0_deriv)
/ (6 * (xInArray[iInterval + 1] - xInArray[iInterval]));
koeffArray[iInterval, 1] = (xInArray[iInterval + 1] * f2x0_deriv - xInArray[iInterval]
* f2x1_deriv) / (2 * (xInArray[iInterval + 1] - xInArray[iInterval]));
koeffArray[iInterval, 2] = ((yInArray[iInterval + 1] - yInArray[iInterval])
- koeffArray[iInterval, 1] * (Math.Pow(xInArray[iInterval + 1], 2)
- Math.Pow(xInArray[iInterval], 2))
- koeffArray[iInterval, 0] * (Math.Pow(xInArray[iInterval + 1], 3)
- Math.Pow(xInArray[iInterval], 3)))
/ (xInArray[iInterval + 1] - xInArray[iInterval]);
koeffArray[iInterval, 3] = yInArray[iInterval] - koeffArray[iInterval, 2] * xInArray[iInterval]
- koeffArray[iInterval, 1] * Math.Pow(xInArray[iInterval], 2)
- koeffArray[iInterval, 0] * Math.Pow(xInArray[iInterval], 3);
}
int start = Convert.ToInt32(xInArray[0]);
int intervalNum = 0;
for (int i = start; i < pointCount + start; i++)
{
while ((Convert.ToDouble(i) >= xInArray[intervalNum] && Convert.ToDouble(i) <= xInArray[intervalNum + 1]) != true)
{
intervalNum++;
}
yOutArray[i - start] = koeffArray[intervalNum, 3] + koeffArray[intervalNum, 2] * i
+ koeffArray[intervalNum, 1] * Math.Pow(i, 2)
+ koeffArray[intervalNum, 0] * Math.Pow(i, 3);
}
}
catch (Exception ex)
{
errorMessage = ex.Message;
return false;
}
return true;
}
