%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Function for calculation/visualisation of an % easy spline example for the 4th exercise. % % A cubic spline for n+1 nodes is calculated % for an equidistant mesh on [x0;xn]. As boundary % conditions, the tangents y'0 and y'n are directly % given (initially zero). % % author: Tobias Neckel % date: 16.11.05 % % INPUT: n number of subintervals in [x0;xn] % setBC string setting the boundary conditions for y'0, y'n % ('zero'==0 or 'deri'==derivative of function) % % OUTPUT: - % % ATTENTION: matlab starts arrays with index 1 instead of 0, % so we have to shift indices when using xi, yi! % function spline_easyN(n, setBC) %close all; %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % calculation numberOfPlotPoints=100; % number of points for eval in subintervals % initialise variables y_prime = zeros(n+1,1); rhs = zeros(n-1,1); % global interval x0=0; xn=4*pi; % mesh size h=(xn-x0)/n; % nodes of mesh xi=(x0 : h : xn)'; yi=eval_f(xi); % set boundary conditions: fix derivatives at x0, xn if setBC=='zero' y_prime(1) = 0; y_prime(end) = 0; tmp_help='0'; elseif setBC=='deri' y_prime(1) = cos(x0); % real derivative of f: cos(x0) y_prime(end) = cos(xn); % real derivative of f: cos(xn) tmp_help='cos(x)'; else error('in spline_easyN: undefined input string for BC!'); end % tridiagonal matrix for linear system of equations % of spline coefficents y'i, i=1,...,n-1 Matrix = 1/h*(diag(4*ones(n-1,1)) + diag(ones(n-2,1),1) + diag(ones(n-2,1), -1) ); % right hand side vector %rhs = 3/h^2*(yi(3:n+1)-yi(1:n-1)); TODO: rauswerfen! % right hand side vector (interior nodes) rhs(2:n-2) = 3/h^2*(yi(4:n)-yi(2:n-2)); % right hand side vector (boundary nodes) rhs(1) = 3/h^2*(yi(3)-yi(1)) - y_prime(1)/h; rhs(n-1) = 3/h^2*(yi(n+1)-yi(n-1)) - y_prime(n+1)/h; % calculate unknown spline coefficients y'i, i=1,...,n-1 y_prime(2:n) = Matrix\rhs; % evaluate spline piecewise t=linspace(0,1,numberOfPlotPoints)'; x=[]; s=[]; [Hzero, Hone, Htwo, Hthree] = HermiteBasis(t); % loop over subintervals i=1,...,n % transformation is here inverse: x calculated by t on % each subinterval for i=1:n %formula: s = y0.*H0(t) + y1.*H1(t) + h*yprime0.*H2(t) + h*yprime1.*H3(t); s = [s; yi(i).*Hzero + yi(i+1).*Hone + h*y_prime(i).*Htwo + h*y_prime(i+1).*Hthree]; x = [x; t*h + xi(i)]; end % calculate maximum pointwise difference error_max = max(abs(s-eval_f(x))); %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % visualisation figure_handle = figure; %%%%%%%%%%%%%%%%%% % plot for matlab plot_handle = plot(x,eval_f(x),'k',x,s,'r',xi,yi,'bo'); set(plot_handle, 'LineWidth', 2); set(gca,'LineWidth',2,'FontWeight','bold'); axis auto; grid on; % zoom on; legend('function f','cubic spline','interpolated nodes',['max. error = ' num2str(error_max)]); %%%%%%%%%%%%%%%%%% % plot for octave %plot(x,eval_f(x),'k;function f;',x,s,'r;cubic spline;',xi,yi,'bo;interpolated nodes;'); title(['Cubic spline through ',num2str(n+1),' interpolation nodes on [0;4\pi] with y_{prime}0 = y_{prime}n = ',tmp_help]); return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Subroutine for evaluation of the function f %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function y=eval_f(x) y=sin(x); return %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Subroutine for evaluation of the Hermite basis % polynomials H0,...,H3 (vector values possible) %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function [Hzero, Hone, Htwo, Hthree]=HermiteBasis(t) Hzero = (1 + 2.*t) .* (1 - t).^2; Hone = (3 - 2.*t) .* (t.^2); Htwo = t .* (1 - t).^2; Hthree = (t - 1) .* (t.^2); return