Autor Tema: Implementación de Splines periódicos

0 Usuarios y 1 Visitante están viendo este tema.

04 Noviembre, 2020, 08:49 pm
Leído 98 veces

FerOliMenNewton

  • $$\Large \color{red}\pi\,\pi\,\pi\,\pi$$
  • Mensajes: 252
  • País: mx
  • Karma: +0/-0
  • Sexo: Masculino
Hola a todos, estoy interesado en implementar en Octave una función que implemente los splines periódicos, como entrada la función tiene vectores \( x,y \) , donde \( x \) son las abcisas de los puntos a interpolar y \( y \) son las ordenadas de los puntos a investigar. La salida es la matriz de coeficientes de los splines \( s_{i}(x)=a_{i}+b_{i}(x-x_{i})+c_{i}(x-x_{i})^2+d_{i}(x-x_{i})^3 \). 
La diferencia con los splines naturales, o los de Malcom-Moller-Forsythe es que estamos interpolando una función \( P \)-periódica en algún intervalo de longitud \( P \), por lo que \( x_{n}-x_{1}=P \) y \( y(1)=y(x_{n}) \), no se exige que las segundas derivadas en los extremos sean cero. Aquí abajo les adjunto el script de Octave que hice, al parecer tengo un error ya que cuando compruebo la matriz de mi programa con la que propone el software estadístico R para los vectores \( x=[-3,-1,0,1,3] \), \(  y=[0.1,0.5,1,0.5,0.1] \) aunque los coeficientes son medianamente apróximados, el error sigue siendo considerable si nos ponemos estricos. Pero ya lo he checado y sigo sin ver mi error.
¿Alguien podría explicarme por qué está mal por favor?
La página en la que estudié el método es la siguiente:
http://nikolavitas.blogspot.com/2013/09/cubic-spline-interpolation-periodic.html
La notación que usa el autor se establece aquí:
http://nikolavitas.blogspot.com/2013/09/cubic-spline-interpolation-natural.html
De antemano muchas gracias :) .
Saludos.

function M=periodicsplines(x,y)
% Implementacion de los splines periodicos:
% x: abscisas de los puntos a interpolar
% y: ordenadas de los puntos a interpolar
% -> matriz A=[a b c] de coeficientes.
n=length(x);
% OBS: Por periodicidad, y debe ser un vector tal que y(1)=y(n).
for i=1:n-1
  h(i)=x(i+1)-x(i);
endfor
for i=1:n-1
  a(i)=y(i);
endfor   
% Construccion de la matriz A del sistema A*x=B , en este caso A sera de tamaño (n-1)*(n-1)
% ya que x=(c1,...,cn-1).
A=zeros(n-1,n-1);
A(1,1)=2*(h(n-1)+h(1));   %2*h(1)
A(1,2)=2*h(1); %2*(h(n-1)+h(1))
A(1,n-1)=h(n-1);
A(n-1,1)=h(n-1);
A(n-1,n-2)=h(n-2);
A(n-1,n-1)=2*(h(n-2)+h(n-1));
for i=2:n-2
A(i,i-1)=h(i-1);
A(i,i)=2*(h(i-1)+h(i));
A(i,i+1)=h(i);
endfor   
% Construccion del vector B
B(1)=(3/h(1))*(a(2)-a(1))-(3/h(n-1))*(a(1)-a(n-1));
B(n-1)=(3/h(n-1))*(a(1)-a(n-1))-(3/h(n-2))*(a(n-1)-a(n-2));
for i=2:n-2
  B(i)=(3/h(i))*(a(i+1)-a(i))-(3/h(i-1))*(a(i)-a(i-1));
endfor
% Resolucion del sistema:
c=inv(A)*B';
% Determinacion de los coeficientes restantes:
a(n)=a(1);
c(n)=c(1);
for j=1:n-1
  b(j)=(1/h(j))*(a(j+1)-a(j))-(h(j)/3)*(2*c(j)+c(j+1));
  d(j)=(c(j+1)-c(j))/(3*h(j));
endfor
% Matriz de coeficientes:
M=[ a(1:n-1)' b' c(1:n-1) d'];
endfunction



05 Noviembre, 2020, 02:40 am
Respuesta #1

FerOliMenNewton

  • $$\Large \color{red}\pi\,\pi\,\pi\,\pi$$
  • Mensajes: 252
  • País: mx
  • Karma: +0/-0
  • Sexo: Masculino
Hola de nuevo,
Ya encontré el error  :P , era un error muy tonto, si se fijan la entrada \( A(1,2) \) es incorrecta, el \( 2 \) no debe ir ahí, no sé como es que no lo vi durante varias horas! Pero venga, así pasa.
Dejaré el código aquí en caso de que le sirva a alguien en el futuro.
Saludos.