% Function implementing natural spline (follow algorith 3.4) function natural_spline %nodes=[0.5, 1, 2, 3, 4, 8,10]; %nodes= [0:.2:10]; nodes=[0.1,0.2,0.3,0.4]; %fnodes=f(nodes); fnodes=[-0.62049958, -0.28398668, 0.00660095, 0.24842440]; %%%%%%%%%%%%%%%%%%%%%%% % Calculate Coefficients of the natural spline % given vaue of the function (fnodes) and the nodes %%%%%%%%%%%%%%%%%%%%%%% n=size(nodes); n=n(2); %determine number of the nodes h = nodes(2:n)-nodes(1:n-1); %create space for meshsteps : h(i)=nodes(i+1)-nodes(i); a=fnodes; % calculate coefficients a_{i} (coinside with the value of the function at nodes) c=zeros(1, n); % setup space for coefficients c b=zeros(1, n-1); % setup space for coefficients b d=zeros(1,n-1); % setup space for coefficients d %%%%%%%%%%%%%%%%%%%%%% % calculate coefficients c_{i}, b_{i}, d_{i} %%%%%%%%%%%%%%%%%%%%%% alpha=(3./h(2:n-1)).*(a(3:n)-a(2:n-1))-(3./h(1:n-2)).*(a(2:n-1)-a(1:n-2)); %% calculate the right side of the system of equations for coefficients c_{i} %% Solving a tridiagonal linear system by algorithm 6.7 l=zeros(1,n); l(1)=1; mu=zeros(1,n-1); z=zeros(1,n-1); % set up coefficients l, mu, z (step 3) for i=2:1:n-1 l(i)=2*(nodes(i+1)-nodes(i-1))-h(i-1)*mu(i-1); mu(i)=h(i)/l(i); z(i)=(alpha(i-1)-h(i-1)*z(i-1))/l(i); end %% Step 5: l(n)=1; c(n)=0; %z(n)=0 %% Step 6: for j=n-1:-1:1 c(j)=z(j)-mu(j)*c(j+1); end for j=1:1:n-1 b(j)=(a(j+1)-a(j))/h(j)-h(j)*(c(j+1)+2*c(j))/3; d(j)=(c(j+1)-c(j))/(3*h(j)); end spline(nodes,a,b,c,d,0.25) spline_der(nodes,b,c,d,0.25) %%%%%%%%%%%%%%%%%%%%%%%%% % Plot the original function and the spline %%%%%%%%%%%%%%%%%%%%%%% figure(1); subplot(1,1,1); clf; % clear the plot % Plot the exact function x=[0:.01:10]; %plot(x, f(x), '-r'); %hold on; % to keep the plot on % calculate the spline m=size(x); m=m(2); y=zeros(1,m); for i=1:1:m y(i)=spline(nodes,a,b,c,d,x(i)); end plot(x,y-f(x),'-b'); hold off; %%%%%%%%%%%%%%%%%%%%%%%%% % at each interval [x_{j} ,x_{j+1}] spline is calculated by the formula % a(j) + b(j)(x-x(j)) + c(j)(x-x(j))^2 + d(j)(x-x(j))^3 % is the point where the spline needs to be evaluated % %%%%%%%%%%%%%%%%%%%%%% function y=spline(nodes,a,b,c,d,x); y=0; n=size(nodes); n=n(2); %determine number of the nodes for j=1:1:n-1 if (x>=nodes(j)) & (x<=nodes(j+1)) y=a(j) + b(j)*(x-nodes(j)) + c(j)*(x-nodes(j))^2 + d(j)*(x-nodes(j))^3; end end function y=spline_der(nodes,b,c,d,x); y=0; n=size(nodes); n=n(2); %determine number of the nodes for j=1:1:n-1 if (x>=nodes(j)) & (x<=nodes(j+1)) y= b(j) + 2*c(j)*(x-nodes(j)) + 3*d(j)*(x-nodes(j))^2; end end function y=f(x) %y=log(x); y=sin(x-5);