Não conheço algoritmo de spline amplamente estudado com esse propósito, portanto tenho a intenção de desenvolver um novo aos poucos.

O que me motivou a perguntar, resolver e ainda registrar aqui no stackoverflow foi eu não conseguir encontrar nem criar um algoritmo eficiente sem estudo aprofundado como esse. De fato, precisei disso porque tive problemas de desempenho por trabalhar com problemas muito grandes com enorme complexidade de operações dos algoritmos. Agora o resultado ficará aqui registrado para quando precisarmos revisar.

Modelando o problema

Para começar, já que:

• fórmulas fm isoladas podem atender a M condições conservando continuidade em derivadas de todas as ordens e
• polinômios m<=n na forma descrita no enunciado de ordens estão dentro do modelo fn ao considerar zero os coeficientes c(m+1), c(m+2), ..., cn (de graus m+1, m+2, ..., n),

então com uma só fórmula fm o spline está completo e atende aos requisitos desejados se m<=n (caso não haja o rigor de não aceitar ordem m de polinômio menor que a ordem n exigida, o que considero verdadeiro) e atender às condições de associação de argumentos e valores de spline (fm(x0)=y0, fm(x1)=y1, ..., fm(xm)=ym).

Resolvido nessas circunstâncias, podemos nos focar então em problemas onde n>m e portanto uma só fórmula não atende a todas as condições de associação. Vamos considerar que cada trecho do spline se limita a argumentos entre dois valores vizinhos de x. O algoritmo desejado atende às seguintes condições.

• Algoritmo tem os parâmetros M ou m, N ou n, x[0, 1, ..., m] e y[0, 1, ..., m].
• Tem N=n+1 coeficientes na fórmula f(x) = polinômio de grau n (não necessariamente o fn do enunciado da questão, deve ser possível resolver os sistemas de equações).
• Resultado é uma função condicional (spline) que associa cada um de m intervalos (X0=[x0,x1], X1=[x1,x2], ..., X(m-1)=[x(m-1),xm], portanto F=m-1 fronteiras x1, x2, ..., x(m-1)) a uma fórmula com modelo f e coeficientes c já calculados.
• Spline tem M=m+1 associações de argumentos e valores de spline (x0->y0, x1->y1, ..., xm->ym) respeitados pelas associações de intervalos às fórmulas com coeficientes adequados.
• Spline tem pelo menos até ordem D=n-1 de derivadas contínuas (f, f', f'', ..., f^(D)) mesmo nas fronteiras.
• Portanto considera-se {D,F,m,M,n,N} c {Inteiros}, sequência x crescente, 0 <= D < n < N <= m < M e F < m.

Leve em conta que não necessariamente deverá se restringir os cálculos ao intervalo [x0,xm], portanto mesmo durante a construção do algoritmo o primeiro intervalo sendo [x0,x1] e o último [x(m-1),xm] na prática o primeiro intevalo do spline final é [-inf,x1] e o último é [x(m-1),+inf].

Reduzindo o problema: pontos iniciais e finais de intervalos

Sendo assim, para cada intervalo Xi = [xi,x(i+1)] = [xi,xf] podemos calcular os N=n+1 coeficientes de f específicos do intervalo. Respeitando as condições f(xi)=yi e f(xf) = y(i+1) = yf, encontramos dois em função dos outros, mas restarão outros n-1 coeficientes a serem calculados em cada trecho.

Se por conveniência escrevermos o polinômio numa forma incomum, f(x) = c0*(xf-x) + (x-xi)*( c1 + (xf-x)*( c2 + c3*x + c4*x² + c5*x³ + ... + cn*x^(n-2) ) ) (por mais confuso que pareça), e resolvermos c0 e c1 com essas duas condições significa a solução do sistema linear { f(xi)=yi , f(x(i+1))=y(i+1) } ser simplesmente { c0=yi/(xf-xi) , c1=yf/(xf-xi) }.

Com isso, o problema se resume em a cada intervalo Xi achar esses n-1 coeficientes c2, c3, ..., cn de f(x) = yi*((xf-x)/(xf-xi)) + (x-xi)*( yf/(xf-xi) + (xf-x)*( c2 + c3*x + c4*x² + c5*x³ + ... + cn*x^(n-2) ) ) respeitando o critério das fronteiras de trechos mantiverem a continuidade de derivadas nas fronteiras.

Reduzindo o problema: modelando o problema das fronteiras

Se há m intervalos (com D=n-1 coeficientes cada a serem calculados), então há F=m-1 fronteiras para ajustar continuidades de D=n-1 derivadas, portanto D*m-D*F coeficientes restantes a serem solucionados, ou seja, restarão D coeficientes, o equivalente aos de um intervalo (como se um ficasse de fora).

Precisaremos adaptar o problema a essa realidade. Uma maneira de fazer isso é criar mais um conjunto de equações que façam o papel de uma fronteira adicional para "preencher as condições do intervalo restante". Outra maneira é substituir os D*m coeficientes por outros valores que correspondam a D*F coeficientes.

Além disso, cada intervalo tem fórmula f com coeficientes de valores distintos das outras e precisamos diferenciá-los. Para isso, cada intervalo Xi terá a fórmula fi com coeficientes ci2, ci3, ..., cin.

Decidi reduzir o número de coeficientes. O modo de redução que escolhi foi para cada uma das D=n-1 ordens o=2,3,...,n substituir os coeficientes c por k de uma das duas seguintes maneiras.

• { c0o=k0o+k1o , c1o=k1o+k2o , ... , c(m-1)o=k(m-1)o+kmo }.
• { c0o=k0o-k1o , c1o=k1p-k2o , ... , c(m-1)o=k(m-1)o-kmo }.

Sendo que em ambas opções considera-se k0o = kmo = 0, restando k1o, k2o, ..., k(m-1)o indefinidos, resultando em F coeficientes por ordem o e assim totalizando F*D. Sendo assim, as fórmulas fi seguem um dos modelos seguintes.

• fi(x) = yi*(x(i+1)-x)/(x(i+1)-xi) + (x-xi)*( y(i+1)/(x(i+1)-xi) + (x(i+1)-x)*( (ki2+k(i+1)2) + (ki3+k(i+1)3)*x + (ki4+k(i+1)4)*x² + (ki5+k(i+1)5)*x³ + ... + (kin+k(i+1)n)*x^(n-2) ) ).
• fi(x) = yi*((x(i+1)-x)/(x(i+1)-xi)) + (x-xi)*( y(i+1)/(x(i+1)-xi) + (x(i+1)-x)*( (ki2-k(i+1)2) + (ki3-k(i+1)3)*x + (ki4-k(i+1)4)*x² + (ki5-k(i+1)5)*x³ + ... + (kin-k(i+1)n)*x^(n-2) ) ).

Com isso, são F=m-1 fronteiras x1, x2, ..., xF e F*D=(m-1)*(n-1) coeficientes k. Para cada fronteira, criar D equações para igualar as derivadas de ordem 1, 2, ..., D leva a um sistema de F*D equações lineares, correspondendo ao número de coeficientes.

Para facilitar, as fórmulas foram implementadas no maple usando o número M de associações de argumentos e resultados, a ordem polinomial n, o índice i do intervalo e o fator kf = ±1 do segundo coeficiente k que compõe coeficiente c (já zerando coeficientes k nulos dos extremos).

Por exemplo, quando M=5 e n=3 as duas fórmulas de fi para cada um dos n=4 intervalos (i=0,1,2,3) são as que estão logo abaixo.

Reduzindo o problema: montando o sistema

Sendo assim, para cada fronteira xi e { x1, x2, ..., xF } tem-se as equações do sistema f(i-1)'(xi) = fi'(xi), f(i-1)''(xi) = fi''(xi), ..., f(i-1)^(D)(xi) = fi^(D)(xi). Para facilitar, foi implementado um procedimento que forma para cada fronteira as equações já com as variáveis e coeficientes reunidos no lado esquedo da equação e os termos independentes a direita, ou seja, o sistema de equações lineares da forma mais legível que pode (ainda que nem tanto quando o problema resolvido é grande) e exibe "fronteira, ordem de derivada, equação correspondente". Lembrando que D = n-1 e F = m-1 = M-2.

Por exemplo, os dois sistemas no exemplo anterior com M=5 e n=3 são exibidos da forma a seguir.

Reduzindo o problema: entendendo o sistema

Alterando o procedimento, reparamos que a representação dos sistemas em forma de matriz parece ser de uma matriz esparsa (preenchida com muitos zeros) segundo algum padrão. No exemplo, as matrizes ficam assim se as colunas se associarem às variáveis na sequência ki_o e { k1_2, k1_3, k2_2, k2_3, k3_2, k3_3 }.

Verificando se cada coeficiente é zero ou outra expressão, podemos observar a esparsidade de matrizes e isso é mais facilmente feito se for exibido na célula ao invés da expressão algum marcador curto que indica se é zero ou não. Alterando o procedimento de modo que isso seja aplicado, observamos algo parecido com o que é apresentado a seguir.

Assim pode-se observar matrizes bem maiores e possivelmente entender o padrão. Além disso, lembremos que cada equação iguala ou subtrai duas fórmulas (uma do intervalo a esquerda da fronteira e a outra da direita), cada uma com até 2*D coeficientes e, ainda assim, metade deles são em comum com a outra, totalizando até 3*D coeficientes k tratadas como variáveis na equação do sistema.

Como será então a matriz nessas condições? Como será quando M=9 e n=2? Será que encontramos a cada linha no máximo três (3*D = 3*(n-1) = 3) céluas não nulas e próximas?

Hmm... aparentemente, sim. E se M=9 e n=3? No máximo seis células não nulas por linha?

E se M=9 e n=4?

Aparentemente, as matrizes somente têm células a cada linha que sejam não nulas numa distância menor que 3*D células a partir da que está na diagonal principal. Pode-se aplicar algoritmos genéricos como o uso de matrizes inversas ou obtenção de forma escalonada reduzida, mas a complexidade de operações é cúbica em relação ao parâmetro F e adaptar um algoritmo que evita trabalhar com todas as células diminui a complexidade em relação a F, podendo assim trabalhar muito melhor com grande número de associação de argumentos e valores de spline principalmente em baixas ordens de polinômios.

Resolvendo o problema

Percebe-se que há F "blocos" de D linhas onde por fora não é esparso. A razão disso é o fato de cada bloco corresponder com uma fronteira, que tem D linhas e trabalha com 2*D ou 3*D coeficientes. O primeiro e o último bloco têm 2*D colunas e os outros têm 3*D.

Aparentemente, é possível separar os blocos como se fossem sistemas distintos. Seria como resolver os sistemas encontrando os valores das variáveis associadas às diagonais principais da matriz inteira, cada uma em função das outras variáveis que se incluem no sistema do bloco. Isso significa 2 sistemas de D equações e 2*D variáveis e F-2 sistemas de D equações e 3*D variáveis, resultando na complexidade F*D³.

Além disso, as D variáveis obtidas de um sistema podem ser usados no próximo para reduzir D colunas (operações de complexidade D² ocorrendo em todos os blocos diferentes do primeiro, logo complexidade F*D²), significando F-1 sistemas de D equações e 2*D variáveis e 1 sistema (o último) de D equações e 3*D variáveis, conservando a complexidade F*D³ mas ainda assim melhorando o desempenho.

Como o último sistema é DxD, espera-se que se obtenha os valores das últimas variáveis. Como os outros são Dx(2*D) e têm valores das variáveis do próximo, é possível obter um sistema DxD delas com operações D² por bloco de maneira similar ao passo anterior.

Isso corresponde com o escalonamento da matriz para a forma escalonada reduzida da matriz onde as operações de linha se restringem às linhas do mesmo bloco e do próximo na primeira etapa e às linhas do mesmo bloco e do anterior na segunda etapa. Essa simples adaptação resolve e reduz a complexidade de (F*D)³ para F*D³.

Com o sistema resolvido, tem-se os coeficientes k, com isso tem-se coeficientes c e finalmente tem-se cada fórmula de cada intervalo com coeficientes calculados, assim obtendo a spline em complexidade F*D³ atendendo às condições de continuidade e associação de argumentos e valores de spline.

Implementando e testando

Usando o maple, pude implementar uma função completa de spline usando kf=1 e testar plotando o gráfico de spline (e suas derivadas) construído com n=3, x=1,3,5,7,9 e y=1,1,2,4,6. O nome escolhido do argumento do spline foi arg. Atente-se ao fato dos nomes de incógnitas estarem diferentes e a abordagem de solução de sistemas escolhida foi simbólica e não matricial.

Observa-se não só neste exemplo mas em outros testes que fiz que os valores estão corretos e as derivadas que devem ser contínuas de fato são. O resultado foi satisfatório. Além disso, fiz uma modificação no algoritmo (passei a tratar k0o = z e usei esse parâmetro z para suavizar a curva) e refiz os mesmos testes. Veja o exemplo anterior com o novo algoritmo.

A suavização ocorre encontrando o z que equilibra os valores absolutos das derivadas de ordem n do primeiro intervalo e do último. Como resultado, encontramos uma nova curva que parece até mais realista que a anterior. Enquanto não encontro problemas com este algoritmo, considero-o satisfatório.