Cubic spline program

Question

I am trying to write a cubic spline interpolation program. I wrote a program, but the schedule does not work out correctly. The spline uses natural boundary conditions (the second dervative at the beginning / end of the node is 0). The code is in Matlab and shown below.

clear all %Function to Interpolate k = 10; %Number of Support Nodes-1 xs(1) = -1; for j = 1:k xs(j+1) = -1 +2*j/k; %Support Nodes(Equidistant) end; fs = 1./(25.*xs.^2+1); %Support Ordinates x = [-0.99:2/(2*k):0.99]; %Places to Evaluate Function fx = 1./(25.*x.^2+1); %Function Evaluated at x %Cubic Spline Code(Coefficients to Calculate 2nd Derivatives) f(1) = 2*(xs(3)-xs(1)); g(1) = xs(3)-xs(2); r(1) = (6/(xs(3)-xs(2)))*(fs(3)-fs(2)) + (6/(xs(2)-xs(1)))*(fs(1)-fs(2)); e(1) = 0; for i = 2:k-2 e(i) = xs(i+1)-xs(i); f(i) = 2*(xs(i+2)-xs(i)); g(i) = xs(i+2)-xs(i+1); r(i) = (6/(xs(i+2)-xs(i+1)))*(fs(i+2)-fs(i+1)) + ... (6/(xs(i+1)-xs(i)))*(fs(i)-fs(i+1)); end e(k-1) = xs(k)-xs(k-1); f(k-1) = 2*(xs(k+1)-xs(k-1)); r(k-1) = (6/(xs(k+1)-xs(k)))*(fs(k+1)-fs(k)) + ... (6/(xs(k)-xs(k-1)))*(fs(k-1)-fs(k)); %Tridiagonal System i = 1; A = zeros(k-1,k-1); while i < size(A)+1; A(i,i) = f(i); if i < size(A); A(i,i+1) = g(i); A(i+1,i) = e(i); end i = i+1; end for i = 2:k-1 %Decomposition e(i) = e(i)/f(i-1); f(i) = f(i)-e(i)*g(i-1); end for i = 2:k-1 %Forward Substitution r(i) = r(i)-e(i)*r(i-1); end xn(k-1)= r(k-1)/f(k-1); for i = k-2:-1:1 %Back Substitution xn(i) = (r(i)-g(i)*xn(i+1))/f(i); end %Interpolation if (max(xs) <= max(x)) error('Outside Range'); end if (min(xs) >= min(x)) error('Outside Range'); end P = zeros(size(length(x),length(x))); i = 1; for Counter = 1:length(x) for j = 1:k-1 a(j) = x(Counter)- xs(j); end i = find(a == min(a(a>=0))); if i == 1 c1 = 0; c2 = xn(1)/6/(xs(2)-xs(1)); c3 = fs(1)/(xs(2)-xs(1)); c4 = fs(2)/(xs(2)-xs(1))-xn(1)*(xs(2)-xs(1))/6; t1 = c1*(xs(2)-x(Counter))^3; t2 = c2*(x(Counter)-xs(1))^3; t3 = c3*(xs(2)-x(Counter)); t4 = c4*(x(Counter)-xs(1)); P(Counter) = t1 +t2 +t3 +t4; else if i < k-1 c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1)); c2 = xn(i+1)/6/(xs(i+1)-xs(i-1+1)); c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6; c4 = fs(i+1)/(xs(i+1)-xs(i-1+1))-xn(i+1)*(xs(i+1)-xs(i-1+1))/6; t1 = c1*(xs(i+1)-x(Counter))^3; t2 = c2*(x(Counter)-xs(i-1+1))^3; t3 = c3*(xs(i+1)-x(Counter)); t4 = c4*(x(Counter)-xs(i-1+1)); P(Counter) = t1 +t2 +t3 +t4; else c1 = xn(i-1+1)/6/(xs(i+1)-xs(i-1+1)); c2 = 0; c3 = fs(i-1+1)/(xs(i+1)-xs(i-1+1))-xn(i-1+1)*(xs(i+1)-xs(i-1+1))/6; c4 = fs(i+1)/(xs(i+1)-xs(i-1+1)); t1 = c1*(xs(i+1)-x(Counter))^3; t2 = c2*(x(Counter)-xs(i-1+1))^3; t3 = c3*(xs(i+1)-x(Counter)); t4 = c4*(x(Counter)-xs(i-1+1)); P(Counter) = t1 +t2 +t3 +t4; end end end P = P'; P(length(x)) = NaN; plot(x,P,x,fx) 

When I run the code, the interpolation function is not symmetrical and does not converge correctly. Can anyone suggest any suggestions about problems in my code? Thanks.

Answer 1

I wrote a cubic spline package in Mathematica a long time ago. Here is my translation of this package in Matlab. Note. I did not consider cubic splines after about 7 years, so I base this on my own documentation. You must check everything that I say.

The main problem is to give us n data points (x(1), y(1)) , ... , (x(n), y(n)) , and we want to calculate piecewise cubic interpolation. Interpolation is defined as

  S(x) = { Sk(x) when x(k) <= x <= x(k+1) { 0 otherwise 

Here Sk (x) is a cubic polynomial of the form

  Sk(x) = sk0 + sk1*(xx(k)) + sk2*(xx(k))^2 + sk3*(xx(k))^3 

Spline Properties:

• The spline passes through the data point Sk(x(k)) = y(k)
• The spline is continuous at the end points and, therefore, continuous everywhere in the interpolation interval Sk(x(k+1)) = Sk+1(x(k+1))
• The spline has a continuous first derivative Sk'(x(k+1)) = Sk+1'(x(k+1))
• The spline has a continuous second derivative Sk''(x(k+1)) = Sk+1''(x(k+1))

To construct a cubic spline from a set of data points, we need to solve for the coefficients sk0 , sk1 , sk2 and sk3 for each of the cubic polynomials n-1 . This is only 4*(n-1) = 4*n - 4 unknowns. Property 1 provides the restrictions n , and properties 2,3,4 each provide additional restrictions n-2 . Thus, we have n + 3*(n-2) = 4*n - 6 restrictions and 4*n - 4 unknowns. This leaves two degrees of freedom. We fix these degrees of freedom, setting the second derivative equal to zero at the initial and final nodes.

Let m(k) = Sk''(x(k)) , h(k) = x(k+1) - x(k) and d(k) = (y(k+1) - y(k))/h(k) . The following is a three-term recurrence relation

  h(k-1)*m(k-1) + 2*(h(k-1) + h(k))*m(k) + h(k)*m(k+1) = 6*(d(k) - d(k-1)) 

m (k) are unknown for which we want to solve. The parameters h(k) and d(k) determined by the input data. This three-term recurrence relation defines a three-diagonal linear system. After determining m(k) coefficients of Sk determined by the expression

  sk0 = y(k) sk1 = d(k) - h(k)*(2*m(k) + m(k-1))/6 sk2 = m(k)/2 sk3 = m(k+1) - m(k)/(6*h(k)) 

Well, this is all the math you need to know in order to fully define the cubic spline calculation algorithm. Here it is in Matlab:

 function [s0,s1,s2,s3]=cubic_spline(x,y) if any(size(x) ~= size(y)) || size(x,2) ~= 1 error('inputs x and y must be column vectors of equal length'); end n = length(x) h = x(2:n) - x(1:n-1); d = (y(2:n) - y(1:n-1))./h; lower = h(1:end-1); main = 2*(h(1:end-1) + h(2:end)); upper = h(2:end); T = spdiags([lower main upper], [-1 0 1], n-2, n-2); rhs = 6*(d(2:end)-d(1:end-1)); m = T\rhs; % Use natural boundary conditions where second derivative % is zero at the endpoints m = [ 0; m; 0]; s0 = y; s1 = d - h.*(2*m(1:end-1) + m(2:end))/6; s2 = m/2; s3 =(m(2:end)-m(1:end-1))./(6*h); 

Here is some code to build a cubic spline:

 function plot_cubic_spline(x,s0,s1,s2,s3) n = length(x); inner_points = 20; for i=1:n-1 xx = linspace(x(i),x(i+1),inner_points); xi = repmat(x(i),1,inner_points); yy = s0(i) + s1(i)*(xx-xi) + ... s2(i)*(xx-xi).^2 + s3(i)*(xx - xi).^3; plot(xx,yy,'b') plot(x(i),0,'r'); end 

Here is a function that builds a cubic spline and enters the well-known Runge function:

 function cubic_driver(num_points) runge = @(x) 1./(1+ 25*x.^2); x = linspace(-1,1,num_points); y = runge(x); [s0,s1,s2,s3] = cubic_spline(x',y'); plot_points = 1000; xx = linspace(-1,1,plot_points); yy = runge(xx); plot(xx,yy,'g'); hold on; plot_cubic_spline(x,s0,s1,s2,s3); 

You can see it in action by doing the following at the Matlab prompt

  >> cubic_driver(5) >> clf >> cubic_driver(10) >> clf >> cubic_driver(20) 

By the time you have twenty nodes, your interpolator is visually indistinguishable from the Runge function.

Some comments on Matlab code: I do not use any for or while loops. I can vectorize all operations. I quickly form a sparse tridiagonal matrix with spdiags . I solve it with the backslash operator. I count on Tim Davis UMFPACK to deal with decomposition, and solve it back and forth.

Hope this helps. The code is available as gistub https://gist.github.com/1269709

Answer 2

An error was detected in the spline function generated by (n-2) to (n-2) should be symmetrical:

 lower = h(2:end); main = 2*(h(1:end-1) + h(2:end)); upper = h(1:end-1); 

http://www.mpi-hd.mpg.de/astrophysik/HEA/internal/Numerical_Recipes/f3-3.pdf