//
// splines cubiques pour une fonnction usuelle
//
function s=reconst(t)
//
// reconstruction de la fonction spline sur [x(i),x(i+1)]
//
i=0;
while t>x(i+1)
i=i+1
end
if (t==xmin) then
   s=f(x(1))
elseif (t==xmax)
   s=f(x(n+2));
else
  xi=t-x(i);
  s=f(x(i))*(2*xi^3-3*h*xi^2+h^3)/h^3+f(x(i+1))*(-2*xi^3+3*h*xi^2)/h^3+...
    p(i)*(xi^3-2*h*xi^2+h^2*xi)/h^2+p(i+1)*(xi^3-h*xi^2)/h^2;	
end
endfunction
//
//  fonction a interpoler
//
xmin=-1;xmax=1;deff('y=f(x)','y=1/(1+x^2)');
n=evstr(x_dialog('nombre de points d''interpolation?','6'));
//
//points d'interpolation
//
n=n-2;
h=(xmax-xmin)/(n+1);
x=xmin:h:xmax;
//
//matrice du systeme a resoudre
//
v=zeros(1,n-2);
A=toeplitz([4,1,v]);
//
// second membre du systeme
//
b=zeros(n,1);
c0=(f(x(2))-f(x(1)))/h;
c1=(f(x(n+2))-f(x(n+1)))/h;
b(1)=(3/h)*(f(x(3))-f(x(1)))-c0;
for i=2:n-1
	b(i)=(3/h)*(f(x(i+2))-f(x(i)));
end;
b(n)=(3/h)*(f(x(n+2))-f(x(n)))-c1;
//
//resolution 1 (routine Scilab)
//
//p=A\b;
// 
//
//resolution 2 (gradient à pas constant)
//
p=ones(n,1);rho=1/6;
for i=1:200
    p=p-rho*(A*p-b)
end
//
p=[c0;p;c1];
//
// reconstruction du profil
//
t=[];sp=[];ft=[];
for i=1:101
  t(i)=xmin+(xmax-xmin)*(i-1)/(100);
  sp(i)=reconst(t(i));
  ft(i)=f(t(i));
end
//
// trace
//
clf()
plot2d([t,t],[ft,sp],leg='courbe exact@courbe interpolee par splines');
//

