// Numerical Algorithms for Visual Computing I - Summer term 2008
// implementation by Kai Hagenburg - May 2008


//interpolation function

function f = interpolate(c,x,p)
[m,n] = size(x)
if n==1
    n = m
 end
  f = c(n)*(p-x(n-1))+ c(n-1)
  for i=2:n-1
    f = f*(p-x(n-i))+c(n-i)
  end
endfunction


//Divided differences
function d = dividiff(x,y)
[n,m]=size(y)
if n==1
    n = m
 end
 d= zeros(n,n)
 d(:,1)=y'
 for j=2:n
  for i=1:n-j+1
   d(i,j) = (d(i+1,j-1)-d(i,j-1))/(x(i+j-1)-x(i))
  end
 end
endfunction


function l=exercise1(x,y,p)
A=dividiff(x,y)
c=A(1,:)
l=interpolate(c,x,p)
endfunction




function fp = newton_interpolation(x,y,p)

n = length(x);
a(1) = y(1);
for k = 1 : n - 1
   d(k, 1) = (y(k+1) - y(k))/(x(k+1) - x(k));
end
for j = 2 : n - 1
   for k = 1 : n - j
      d(k, j) = (d(k+1, j - 1) - d(k, j - 1))/(x(k+j) - x(k));
   end
end
d
for j = 2 : n
   a(j) = d(1, j-1);
end
Df(1) = 1;
c(1) = a(1);
for j = 2 : n
   Df(j)=(p - x(j-1)) .* Df(j-1);
   c(j) = a(j) .* Df(j);
end
fp=sum(c);
endfunction

//====================================================
// Runge function (exercise 2)
function y = f(x)
y=1.0/(x^2+1.0)
endfunction

// sqrt of abs function (exercise 4)
function y = g(x) 
y=sqrt(abs(x))
endfunction


function x=changer(n)
for i=1:n+1
  y(i) = -5*cos((i-1)/n*%pi)
end

if (modulo(n+1,2)==0)
 n_half = (n+1)/2
 x(1) = y(n_half)
 j=1
 for i=1:n_half
  j=j+1
  x(j)=y(n_half+i)
  j=j+1
  if (i<>n_half)
  x(j)=y(n_half-i)
  end
 end
end

if (modulo(n+1,2)==1)
 n_half=ceil((n+1)/2)
 x(1)=y(n_half)
 j=1
 for i=1:n_half-1
  j=j+1
  x(j)=y(n_half+i)
  j=j+1
  x(j)=y(n_half-i)
 end
end
endfunction
 


function t=exercise2a(n)
timer()
// create interpolation abscissae
for i=1:n+1
  y(i) = -5+10/n*(i-1)
end

//create function values for interpolation points
for i=1:n+1
 z(i)=f(y(i))
end

//compute coefficients via Divided differences
A=dividiff(y',z')

//extract coefficients for the interpolation polynomial
//which are stored in the first row of the 
//Divided Difference matrix
c = A(1,:)

//derive polynomial
x_axis=linspace(-5,5,1001)
for i=1:1001
 interpolated_function(i)=interpolate(c,y',x_axis(i))
 original_function(i) = f(x_axis(i))
end
xset("window",0)
plot2d(x_axis,interpolated_function,style=1,leg="interpolation")
plot2d(x_axis,original_function,style=2,leg="original")
t=timer()
endfunction

function exercise2b(n)
// create interpolation abscissae
//for i=1:n+1
//  y(i) = -5*cos((i-1)/n*%pi)
//end

y=changer(n)

//create function values for interpolation points
for i=1:n+1
 z(i)=f(y(i))
end

//compute coefficients via Divided differences
A=dividiff(y',z')

//extract coefficients for the interpolation polynomial
//which are stored in the first row of the 
//Divided Difference matrix
c = A(1,:)

//derive polynomial
x_axis=linspace(-5,5,101)

for i=1:101
 //
 
 interpolated_function(i)=interpolate(c,y',x_axis(i))
 original_function(i) = f(x_axis(i))
end
xset("window",0)
plot2d(x_axis,interpolated_function,style=1,leg="interpolation")
plot2d(x_axis,original_function,style=2,leg="original")
endfunction




function test
// create interpolation abscissae
y(1) = 0
y(2) = 1
y(3) = 2

z(1) = 8
z(2) = 5
z(3) = 4

//create function values for interpolation points

//compute coefficients via Divided differences
A=dividiff(y',z')

//extract coefficients for the interpolation polynomial
//which are stored in the first row of the 
//Divided Difference matrix
c = A(1,:)

//derive polynomial
x_axis=linspace(1,4,101)
//inter_function=poly(c,'x','coeff')

for i=1:101
 interpolated_function(i)=interpolate(c,y,x_axis(i))
 original_function(i) = x_axis(i)^2-4*x_axis(i)+8
 //original_function(i) = f(x_axis(i))
end
xset("window",0)
plot2d(x_axis,interpolated_function,style=1,leg="interpolation")
plot2d(x_axis,original_function,style=2,leg="original")
endfunction


function c = compute_spline(f,d,x)
[n,m]=size(f)
if n==1
    n = m
 end
c(:,1) = f'
c(:,2) = d'
for i=1:n-1
 c(i,3) = (3*f(i+1)-3*f(i)-2*d(i)*(x(i+1)-x(i))-d(i+1)*(x(i+1)-x(i)))/(x(i+1)-x(i))^2
 c(i,4) = (2*f(i)-2*f(i+1)+d(i)*(x(i+1)-x(i))+d(i+1)*(x(i+1)-x(i)))/(x(i+1)-x(i))^3
 end
 
endfunction



function error = exercise4a(n)
// create interpolation abscissae
//for i=1:n+1
//  y(i) = 1*cos((i-1)/n*%pi)
//end

y=changer(n)

//create function values for interpolation points
for i=1:n+1
 z(i)=g(y(i))
end

//compute coefficients via Divided differences
A=dividiff(y',z')

//extract coefficients for the interpolation polynomial
//which are stored in the first row of the 
//Divided Difference matrix
c = A(1,:)

//derive polynomial
x_axis=linspace(-1,1,101)

for i=1:101
 //
 //interpolated_function(i)=interpolate(c,y,x_axis(i))
 interpolated_function(i) = newton_interpolation(y,z,x_axis(i))
 original_function(i) = g(x_axis(i))
end
xset("window",0)
plot2d(x_axis,interpolated_function,style=1,leg="interpolation")
plot2d(x_axis,original_function,style=2,leg="original")

for i=1:n+1
 inter(i) = interpolate(c,y,y(i))
 end
 
error = 1
//max(abs(inter - y))
endfunction

function [A,b] = build_natural_system(x,f)
[n,m]=size(f)
if n==1
    n = m
 end
A=zeros(n,n)
// Boundaries for A and b
A(1,1) = 2
A(1,2) = 1

A(n,n-1) = 1
A(n,n) = 2

b(1) = 3*(f(2)-f(1))/(x(2)-x(1))
b(n) = 3*(f(n)-f(n-1))/(x(n)-x(n-1))

// standard entries
for i=2:n-1
A(i,i) = 2*(x(i+1)-x(i) + x(i)-x(i-1))
A(i,i-1) = x(i+1)-x(i)
A(i,i+1) = x(i)-x(i-1)
b(i) = 3*((f(i+1)-f(i))*(x(i)-x(i-1))/(x(i+1)-x(i))+ (f(i)-f(i-1))*(x(i+1)-x(i))/(x(i)-x(i-1)))
end
endfunction


function [A,b] = build_complete_system(x,f)
[n,m]=size(f)
if n==1
    n = m
 end
A=zeros(n,n)
// Boundaries for A and b
A(1,1) = 1
A(n,n) = 1

b(1) = (f(2)-f(1))/(x(2)-x(1))
b(n) = (f(n)-f(n-1))/(x(n)-x(n-1))

// standard entries
for i=2:n-1
A(i,i) = 2*(x(i+1)-x(i) + x(i)-x(i-1))
A(i,i-1) = x(i+1)-x(i)
A(i,i+1) = x(i)-x(i-1)
b(i) = 3*((f(i+1)-f(i))*(x(i)-x(i-1))/(x(i+1)-x(i))+ (f(i)-f(i-1))*(x(i+1)-x(i))/(x(i)-x(i-1)))
end
endfunction


function error = exercise4b_natural(n)
for i=1:n+1
  x(i) = -1*cos((i-1)/n*%pi)
end

for i=1:n+1
 y(i)=g(x(i))
end

[A,b] = build_natural_system(x',y')
d=A\b
c = compute_spline(y',d',x')


xset("window",0)
x_axis=linspace(-1,1,1001)
for i=1:1001
 original_function(i) = g(x_axis(i))
end
plot2d(x_axis,original_function,style=2,leg="original")
k=1
for i=1:n
  x_axis_1 = linspace(x(i),x(i+1),101)
  for j=1:101
  interpolated_function(j) = c(i,1)+c(i,2)*(x_axis_1(j)-x(i))+c(i,3)*(x_axis_1(j)-x(i))^2 + c(i,4)*(x_axis_1(j)-x(i))^3
  error_inter(k) = interpolated_function(j)
  error_orig(k) = g(x_axis_1(j))
  k=k+1
  end
  compare(i)=interpolated_function(1)
  compare(i+1) = interpolated_function(101)
  plot2d(x_axis_1,interpolated_function,style=1)
end
   error = max (abs(error_inter-error_orig ))

endfunction

function error = exercise4b_complete(n)
for i=1:n+1
  x(i) = -1*cos((i-1)/n*%pi)
end

for i=1:n+1
 y(i)=g(x(i))
end

[A,b] = build_complete_system(x',y')
d=A\b
c = compute_spline(y',d',x')


xset("window",0)
x_axis=linspace(-1,1,1001)
for i=1:1001
 original_function(i) = g(x_axis(i))
end


plot2d(x_axis,original_function,style=2,leg="original")

k=1
for i=1:n
  x_axis_1 = linspace(x(i),x(i+1),101)
  for j=1:101
  interpolated_function(j) = c(i,1)+c(i,2)*(x_axis_1(j)-x(i))+c(i,3)*(x_axis_1(j)-x(i))^2 + c(i,4)*(x_axis_1(j)-x(i))^3
  error_inter(k) = interpolated_function(j)
  error_orig(k) = g(x_axis_1(j))
  k=k+1
  end
  compare(i)=interpolated_function(1)
  compare(i+1) = interpolated_function(101)
  plot2d(x_axis_1,interpolated_function,style=1)
end
   
error = max (abs(error_inter-error_orig ))
endfunction




function  exercise3b_spline(n)
x(1) = -1
x(2) = 0
x(3) = 1
x(4) = 3

y(1) = 1
y(2) = 3
y(3) = 3
y(4) = 45

[A,b] = build_natural_system(x',y')
d=A\b
c = compute_spline(y',d',x')


xset("window",0)
//x_axis=linspace(-1,1,1001)

plot2d(x,y,style=2,leg="original")

for i=1:3
  x_axis_1 = linspace(x(i),x(i+1),101)
  for j=1:101
  interpolated_function(j) = c(i,1)+c(i,2)*(x_axis_1(j)-x(i))+c(i,3)*(x_axis_1(j)-x(i))^2 + c(i,4)*(x_axis_1(j)-x(i))^3
  end
  plot2d(x_axis_1,interpolated_function,style=1)
end
   

endfunction


function [alpha,errorvec,sq] = errorcalc(n)
for i=2:2:n
 errorvec(i) = exercise4b_natural(i)
end

for i=2:2:n
  sq(i)=sqrt(i+1)*errorvec(i)
end

for i=2:2:n-2
alpha(i) = log( errorvec(i+2)/errorvec(i)) / log( (i+1)/(i+3) )
//alpha(i) = log( errorvec(i+2)/errorvec(i)) / log( (2/(i+2))/(2/(i)) )
end
endfunction