// 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