Pasted

#include <stdio.h>
#include <math.h>
 
// FALTA: 'SPLINE (PRECISA SPLINE?)
 
/*
Funções básicas de interpolação, do tipo linear,
bilinear, polinomial e spline.
Para as interpolações lineares, polinomiais e spline,
o funcionamento é:
 
    y0 = interpol_xxx(n, x, y, x0)
 
onde x e y são vetores de double de tamanho n, que 
descrevem em apenas alguns pontos uma função y = f(x).
O x0 é um double que é o ponto que queremos obter a
interpolação em questão.
 
A interpolação bilinear possui o seguinte funcionamento:
 
    z0 = interpol_bilinear(n, m, x, y, z, x0, y0)
 
onde x é um vetor de tamanho n e y é um vetor de tamanho m,
z é uma matriz n por m que descreve uma função z = f(x,y)
em alguns pontos. Importante notar que o algoritmo realiza
a interpolação primeiro no eixo y e depois no x.
 
Exemplo:
 
    Função y(x) = x^2 + x - 1.
    Definimos x = {1,2,3,4,5} e n = 5.
    Logo obtemos y = {1,5,11,19,29}, aplicando y(x) em um
    loop.
    -- A partir de agora, só temos conhecimento de x e y,
    -- a expressão de y(x) não é mais conhecida.
    Queremos y0 = y(x0) em x0 = 3.5. Claramente não temos
    o valor exato nos vetores, logo necessitamos realizar
    uma interpolação:
    x0 = 3.5
    y0 = interpol_xxx(n, x, y, x0)
    No final temos que y0 = 15, muito próximo do resultado
    real, y(3.5) = 14.75
 
Observações:
 
Um polinômio de grau N-1 pode ser descrito perfeitamente por N pontos, e vice-versa.
 
Para muitos pontos (n > 5), a interpolação polinomial se
torna instável, então não é recomendada.
*/
 
double interpol_linear(int n, double x[n], double y[n], double x0)
{
    int k, i;
    for(i = 1; i < n; i++) // Encontrar índices vizinhos
    {
        if(x0 <= x[i])
        {
            k = i - 1;
            break;
        }
    }
    if(i == n) k = i - 2;
 
    return (y[k+1] - y[k])/(x[k+1] - x[k])*(x0 - x[k]) + y[k];
}
 
double interpol_polinomial(int n, double x[n], double y[n], double x0)
{
    double soma = 0;
    int i, j;
 
    for(i = 0; i < n; i++)
    {
        double produto = 1;
        for(j = 0; j < n; j++)
        {
            if(j != i)
            {
                produto *= (x0 - x[j])/(x[i] - x[j]);
            }
        }
        soma += produto * y[i];
    }
 
    return soma;
}
 
double interpol_bilinear(int n, int m, double x[n], double y[m], double z[n][m], double x0, double y0)
{
    double zn1, zn2;
 
    int i,kx;
    for(i = 1; i < n; i++) // Encontrar índices vizinhos
    {
        if(x0 <= x[i])
        {
            kx = i - 1;
            break;
        }
    }
 
    zn1 = interpol_linear(m, y, z[kx], y0);
    zn2 = interpol_linear(m, y, z[kx+1], y0);
 
    double xpos[2] = {x[kx],x[kx+1]};
    double zpos[2] = {zn1,zn2};
    return interpol_linear(2, xpos, zpos,x0);
}
 
/* ----------- COMEÇO DO EXEMPLO ----------- */
 
double f(double x)
{
    return 0.1*x*x*x + x*x + x - 1;
}
 
int main(void)
{
    double x[4] = {1,2,3,4}, y[4], x0 = 0.9;
 
    int i;
    for(i = 0; i < 4; i++) y[i] = f(x[i]);
 
    printf("x\t=\t%f\nf(x)\t=\t%f\nl(x)\t=\t%f\np(x)\t=\t%f\n", x0, f(x0), interpol_linear(4,x,y,x0), interpol_polinomial(4,x,y,x0));
 
    double zn[2][2] ={{0,1},{2,3}};
    double xn[2] = {0,1};
    double yn[2] = {0,1};
    printf("%f",interpol_bilinear(2,2,xn,yn,zn,0.5,0.5));
 
    // Interpolação polinomial retorna ponto exato pois a
    // função é de 3° grau, e são 4 pontos.
}
 
 
/* ------------ FIM DO EXEMPLO ------------ */