cubic_spline_get_coeffs

cubic_spline_get_coeffs - N=length(X)-1; LY=length(Y);...

Info iconThis preview shows pages 1–2. Sign up to view the full content.

View Full Document Right Arrow Icon
function S=cubic_spline_get_coeffs(X,Y) %% Computes the spline coefficients for a %% natural or clamped cubic spline interpolant %% with modifications made for free end conditions %% Usage: %% for the natural cubic spline: S=cubic_spline_get_coeffs(X,Y) %% for clamped cubic spline: S=cubic_spline_get_coeffs(X,[dx0 Y dxn]) %% Input: %% X: the x coordinates for the knots, sorted in ascending order %% Y: the y coordinates for the knots %% if length(Y) < 2 + length(X), it is assumed that %% free end conditions are used. %% if length(Y) >= 2 + length(X), it is assumed that %% clamped end conditions are used. The first Y value is the %% slope at the first knot, the last Y value is the slope at %% the last knot. The y coordinates of the knots are stored in %% Y(2), . .., Y(length(X)+2). %% Output: %% S: the rows of S are the coefficients, %% the cubic for the interval between X(k) and X(k+1) is %% S(k,1)*(x-X(k))^3 + S(k,2)*(x-X(k))^2 + S(k,3)*(x-X(k)) + S(k,4)
Background image of page 1

Info iconThis preview has intentionally blurred sections. Sign up to view the full version.

View Full DocumentRight Arrow Icon
Background image of page 2
This is the end of the preview. Sign up to access the rest of the document.

Unformatted text preview: N=length(X)-1; LY=length(Y); clamped = ((N+3) &lt;= LY); if (clamped) dx0=Y(1); dxn=Y(N+3); Y = Y(2:N+2); end H=diff(X); D=diff(Y)./H; A=H(2:N-1); B=2*(H(1:N-1)+H(2:N)); C=H(2:N); U=6*diff(D); if (clamped) %% Clamped cubic spline endpoint conditions B(1)=B(1)-H(1)/2; U(1)=U(1)-3*(D(1)-dx0); B(N-1)=B(N-1)-H(N)/2; U(N-1)=U(N-1)-3*(dxn-D(N)); end; %% Solve the tridiagonal system of equations %% elimination phase for k=2:N-1 temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end; %% Solve and back-substitute M(N)=U(N-1)/B(N-1); for k=N-2:-1:1 M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end; if (clamped) %% Clamped cubic spline endpoint conditions M(1)=3*(D(1)-dx0)/H(1)-M(2)/2; M(N+1)= 3*(dxn-D(N))/H(N)-M(N)/2; else %% Free spline endpoint conditions M(1)=0; M(N+1)=0; end; %% The spline coefficients for k=0:N-1 S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1)); S(k+1,2)=M(k+1)/2; S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6; S(k+1,4)=Y(k+1); end;...
View Full Document

Page1 / 2

cubic_spline_get_coeffs - N=length(X)-1; LY=length(Y);...

This preview shows document pages 1 - 2. Sign up to view the full document.

View Full Document Right Arrow Icon
Ask a homework question - tutors are online