Códigos sección 3.
% NEVILLE'S ITERATED INTERPOLATION ALGORITHM 3.1
%
% To evaluate the interpolating polynomial P on the
% (n+1) distinct numbers x(0), ..., x(n) at the number x
% for the function f:
%
% INPUT: numbers x(0),..., x(n) as XX(0),...,XX(N);
% number x; values of f as the first column of Q
% or may be computed if function f is supplied.
%
% OUTPUT: the table Q with P(x) = Q(N+1,N+1).
syms('TRUE', 'FALSE', 'OK', 'FLAG', 'N', 'I', 'XX', 'Q', 'A', 'NAME');
syms('INP', 'X', 'D', 'J', 'OUP','x','s');
TRUE = 1;
FALSE = 0;
fprintf(1,'This is Nevilles Method.\n');
OK = FALSE;
while OK == FALSE
fprintf(1,'Choice of input method:\n');
fprintf(1,'1. Input entry by entry from keyboard\n');
fprintf(1,'2. Input data from a text file\n');
fprintf(1,'3. Generate data using a function F\n');
fprintf(1,'Choose 1, 2, or 3 please\n');
FLAG = input(' ');
if FLAG == 1 | FLAG == 2 | FLAG == 3
OK = TRUE;
end
end
if FLAG == 1
OK = FALSE;
while OK ~= TRUE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
OK = TRUE;
XX = zeros(N+1);
Q = zeros(N+1,N+1);
for I = 0:N
fprintf(1,'Input X(%d) and F(X(%d)) ', I, I);
fprintf(1,'on separate lines.\n');
XX(I+1) = input(' ');
Q(I+1,1) = input(' ');
end
else
fprintf(1,'Number must be a positive integer\n');
end
end
end
if FLAG == 2
fprintf(1,'Has a text file been created with the data in two
columns?\n');
fprintf(1,'Enter Y or N\n');
A = input(' ','s');
if A == 'Y' | A == 'y'
fprintf(1,'Input the file name in the form - ');
fprintf(1,'drive:\\name.ext\n');
fprintf(1,'For example: A:\\DATA.DTA\n');
NAME = input(' ','s');
INP = fopen(NAME,'rt');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input N\n');
N = input(' ');
if N > 0
XX = zeros(N+1);
Q = zeros(N+1,N+1);
for I = 0:N
XX(I+1) = fscanf(INP, '%f',1);
Q(I+1,1) = fscanf(INP, '%f',1);
end
fclose(INP);
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n');
end
end
else
fprintf(1,'Please create the input file in two column ');
fprintf(1,'form with the X values and\n');
fprintf(1,'F(X) values in the corresponding columns.\n');
fprintf(1,'The program will end so the input file can ');
fprintf(1,'be created.\n');
OK = FALSE;
end
end
if FLAG == 3
fprintf(1,'Input the function F(x) in terms of x\n');
fprintf(1,'For example: cos(x)\n');
s = input(' ');
F = inline(s,'x');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
XX = zeros(N+1);
Q = zeros(N+1,N+1);
for I = 0:N
fprintf(1,'Input X(%d)\n', I);
XX(I+1) = input(' ');
Q(I+1,1) = F(XX(I+1));
end
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n');
end
end
end
if OK == TRUE
fprintf(1,'Input point at which the polynomial is to be evaluated\n');
X = input(' ');
end
if OK == TRUE
% STEP 1
D = zeros(N+1);
D(1) = X-XX(1);
for I = 1:N
D(I+1) = X-XX(I+1);
for J = 1:I
Q(I+1,J+1) = (D(I+1)*Q(I,J)-D(I-J+1)*Q(I+1,J))/(D(I+1)-D(I-J+1));
end
end
% STEP 2
fprintf(1,'Select output destination\n');
fprintf(1,'1. Screen\n');
fprintf(1,'2. Text file\n');
fprintf(1,'Enter 1 or 2\n');
FLAG = input(' ');
if FLAG == 2
fprintf(1,'Input the file name in the form - drive:\\name.ext\n');
fprintf(1,'For example: A:\\OUTPUT.DTA\n');
NAME = input(' ','s');
OUP = fopen(NAME,'wt');
else
OUP = 1;
end
fprintf(OUP, 'NEVILLES METHOD\n');
fprintf(OUP, 'Table for P evaluated at X = %12.8f , follows: \n', X);
fprintf(OUP, 'Entries are XX(I), Q(I,0), ..., Q(I,I) ');
fprintf(OUP, 'for each I = 0, ..., N where N = %3d\n\n', N);
for I = 0:N
fprintf(OUP, '%11.8f ', XX(I+1));
for J = 0:I
fprintf(OUP, '%11.8f ', Q(I+1,J+1));
end
fprintf(OUP, '\n');
end
if OUP ~= 1
fclose(OUP);
fprintf(1,'Output file %s created successfully\n',NAME);
end
end
% NEWTONS INTERPOLATORY DIVIDED-DIFFERENCE FORMULA ALGORITHM 3.2
% To obtain the divided-difference coefficients of the
% interpolatory polynomial P on the (n+1) distinct numbers x(0),
% x(1), ..., x(n) for the function f:
% INPUT: numbers x(0), x(1), ..., x(n); values f(x(0)), f(x(1)),
% ..., f(x(n)) as the first column Q(0,0), Q(1,0), ...,
% Q(N,0) of Q, or may be computed if function f is supplied.
% OUTPUT: the numbers Q(0,0), Q(1,1), ..., Q(N,N) where
% P(x) = Q(0,0) + Q(1,1)*(x-x(0)) + Q(2,2)*(x-x(0))*
% (x-x(1)) +... + Q(N,N)*(x-x(0))*(x-x(1))*...*(x-x(N-1)).
syms('OK', 'FLAG', 'N', 'I', 'X', 'Q', 'A', 'NAME', 'INP','OUP', 'J');
syms('s','x');
TRUE = 1;
FALSE = 0;
fprintf(1,'Newtons form of the interpolation polynomial\n');
OK = FALSE;
while OK == FALSE
fprintf(1,'Choice of input method:\n');
fprintf(1,'1. Input entry by entry from keyboard\n');
fprintf(1,'2. Input data from a text file\n');
fprintf(1,'3. Generate data using a function F\n');
fprintf(1,'Choose 1, 2, or 3 please\n');
FLAG = input(' ');
if FLAG == 1 | FLAG == 2 | FLAG == 3
OK = TRUE;
end
end
if FLAG == 1
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
OK = TRUE;
X = zeros(N+1);
Q = zeros(N+1,N+1);
for I = 0:N
fprintf(1,'Input X(%d) and F(X(%d)) ', I, I);
fprintf(1,'on separate lines\n');
X(I+1) = input(' ');
Q(I+1,1) = input(' ');
end
else
fprintf(1,'Number must be a positive integer\n');
end
end
end
if FLAG == 2
fprintf(1,'Has a text file been created with the data in two
columns?\n');
fprintf(1,'Enter Y or N\n');
A = input(' ','s');
if A == 'Y' | A == 'y'
fprintf(1,'Input the file name in the form - ');
fprintf(1,'drive:\\name.ext\n');
fprintf(1,'For example: A:\\DATA.DTA\n');
NAME = input(' ','s');
INP = fopen(NAME,'rt');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
X = zeros(N+1);
Q = zeros(N+1,N+1);
for I = 0:N
X(I+1) = fscanf(INP, '%f',1);
Q(I+1,1) = fscanf(INP, '%f',1);
end
fclose(INP);
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n')
end
end
else
fprintf(1,'Please create the input file in two column ');
fprintf(1,'form with the X values and\n');
fprintf(1,'F(X) values in the corresponding columns.\n');
fprintf(1,'The program will end so the input file can ');
fprintf(1,'be created.\n');
OK = FALSE;
end
end
if FLAG == 3
fprintf(1,'Input the function F(x) in terms of x\n');
fprintf(1,'For example: cos(x)\n');
s = input(' ');
F = inline(s,'x');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
X = zeros(N+1);
Q = zeros(N+1,N+1);
for I = 0:N
fprintf(1,'Input X(%d)\n', I);
X(I+1) = input(' ');
Q(I+1,1) = F(X(I+1));
end
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n');
end
end
end
if OK == TRUE
fprintf(1,'Select output destination\n');
fprintf(1,'1. Screen\n');
fprintf(1,'2. Text file\n');
fprintf(1,'Enter 1 or 2\n');
FLAG = input(' ');
if FLAG == 2
fprintf(1,'Input the file name in the form - drive:\\name.ext\n');
fprintf(1,'For example: A:\\OUTPUT.DTA\n');
NAME = input(' ','s');
OUP = fopen(NAME,'wt');
else
OUP = 1;
end
fprintf(OUP, 'NEWTONS INTERPOLATION POLYNOMIAL\n\n');
% STEP 1
for I = 1:N
for J = 1:I
Q(I+1,J+1) = (Q(I+1,J) - Q(I,J)) / (X(I+1) - X(I-J+1));
end
end
% STEP 2
fprintf(OUP, 'Input data follows:\n');
for I = 0:N
fprintf(OUP,'X(%d) = %12.8f F(X(%d)) = %12.8f\n',I,X(I+1),I,Q(I+1,1));
end
fprintf(OUP, '\nThe coefficients Q(0,0), ..., Q(N,N) are:\n');
for I = 0:N
fprintf(OUP, '%12.8f\n', Q(I+1,I+1));
end
if OUP ~= 1
fclose(OUP);
fprintf(1,'Output file %s created successfully\n',NAME);
end
end
Diferencias divididas con gráfica
%Programa para interpolacion por el metodo de las diferencias divididas
% de newton
clear all
n=input('Numero de datos =');
x=input('x= ');
y=input('y= ');
for i=1:n
M(i,1)=y(i);
end
for j=2:n
for i=1:n-j+1
M(i,j)=(M(i+1,j-1)-M(i,j-1))/(x(i+j-1)-x(i));
end
end
%Programa para efectuar la productoria:
%(t+x(1))(t+x(2))(t+x(3))......(t+x(m))
%El resultado es una matriz de la forma:
%n=1: q(1,1) q(1,2)
%n=2: q(2,1) q(2,2) q(2,3)
%n=3: q(3,1) q(3,2) q(3,3) q(3,4)
%n=4: q(4,1) q(4,2) q(4,3) q(4,4) q(4,5)
%m=input('numero de factores =');
%for i=1:m
%x(i)=input(strcat('entre la abscisa x ',num2str(i), ' ='));
%end
% El algoritmo para llenar la matriz es el siguiente:
m=n-1;
X=-x;
q(1,1)=X(1);
for i=1:n
q(i,i+1)=1;
end
for i=2:n
q(i,1)=X(i)*q(i-1,1);
end
for i=2:n
for j=2:n-1
q(i,j)=X(i)*q(i-1,j)+q(i-1,j-1);
end
end
for i=3:n
q(i,i)=X(i)+q(i-1,i-1);
end
%Programa para expandir el polinomio de Newton, el cual queda en la
forma:
%M(1,1)+M(1,2)*(t+X(1))+M(1,3)*(t+X(1))*(t+X(2))
%+M(1,4)*(t+X(1))*(t+X(2))*(t+X(3))+...
+M(1,n)*(t+X(1))*(t+X(2))*(t+X(3))...(t+X(m))
%En su forma expandida es el polinomio de grado:m
%a(1)+a(2)*x+a(3)*(x^2)+a(4)*(x^3)+....+a(n)*(x^m)
a(1)=M(1,1);
for i=2:n
a(1)=a(1)+M(1,i)*q(i-1,1);
end
for i=2:n
a(i)=0;
end
for i=2:n
for j=i:n
a(i)=a(i)+M(1,j)*q(j-1,i);
end
end
c=fliplr(a);
poly2str(c,'t')
t=x(1):0.01:x(n);
p=polyval(c,t);
plot(x,y,'o');
hold on
plot(t,p);
grid on
% HERMITE INTERPOLATION ALGORITHM 3.3
%
% TO OBTAIN THE COEFFICIENTS OF THE HERMITE INTERPOLATING
% POLYNOMIAL H ON THE (N+1) DISTINCT NUMBERS X(0), ..., X(N)
% FOR THE FUNCTION F:
%
% INPUT: NUMBERS X(0), X(1), ..., X(N); VALUES F(X(0)), F(X(1)),
% ..., F(X(N)) AND F'(X(0)), F'(X(1)), ..., F'(X(N)).
%
% OUTPUT: NUMBERS Q(0,0), Q(1,1), ..., Q(2N + 1,2N + 1) WHERE
%
% H(X) = Q(0,0) + Q(1,1) * ( X - X(0) ) + Q(2,2) *
% ( X - X(0) )**2 + Q(3,3) * ( X - X(0) )**2 *
% ( X - X(1) ) + Q(4,4) * ( X - X(0) )**2 *
% ( X - X(1) )**2 + ... + Q(2N + 1,2N + 1) *
% ( X - X(0) )**2 * ( X - X(1) )**2 * ... *
% ( X - X(N - 1) )**2 * (X - X(N) ).
syms('OK', 'FLAG', 'N', 'I', 'X', 'Q', 'A', 'NAME', 'INP');
syms('Z', 'K', 'J', 'OUP', 'XX', 'S','x','s','I1');
TRUE = 1;
FALSE = 0;
fprintf(1,'This is Hermite interpolation.\n');
OK = FALSE;
while OK == FALSE
fprintf(1,'Choice of input method:\n');
fprintf(1,'1. Input entry by entry from keyboard\n');
fprintf(1,'2. Input data from a text file\n');
fprintf(1,'3. Generate data using a function F\n');
fprintf(1,'Choose 1, 2, or 3 please\n');
FLAG = input(' ');
if FLAG == 1 | FLAG == 2 | FLAG == 3
OK = TRUE;
end
end
if FLAG == 1
OK = FALSE;
while OK == FALSE
fprintf(1,'Input the number of data points minus 1\n');
N = input(' ');
if N > 0
OK = TRUE;
X = zeros(N+1);
Q = zeros(2*N+2,2*N+2);
for I = 0:N
fprintf(1,'Input X(%d), F(X(%d)), and ', I, I);
fprintf(1,'F''(X(%d)) on separate lines\n ', I);
X(I+1) = input(' ');
Q(2*I+1,1) = input(' ');
Q(2*I+2,2) = input(' ');
end
else
fprintf(1,'Number must be a positive integer\n');
end
end
end
if FLAG == 2
fprintf(1,'Has a text file been created with the data in three
columns?\n');
fprintf(1,'Enter Y or N\n');
A = input(' ','s');
if A == 'Y' | A == 'y'
fprintf(1,'Input the file name in the form - ');
fprintf(1,'drive:\\name.ext\n');
fprintf(1,'for example: A:\\DATA.DTA\n');
NAME = input(' ','s');
INP = fopen(NAME,'rt');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input the number of data points minus 1\n');
N = input(' ');
if N > 0
X = zeros(N+1);
Q = zeros(2*N+2,2*N+2);
for I = 0:N
X(I+1) = fscanf(INP, '%f',1);
Q(2*I+1,1) = fscanf(INP, '%f',1);
Q(2*I+2,2) = fscanf(INP, '%f',1);
end
fclose(INP);
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n');
end
end
else
fprintf(1,'Please create the input file in three column ');
fprintf(1,'form with the X values, F(X), and\n');
fprintf(1,'derivative values in the corresponding columns.\n');
fprintf(1,'The program will end so the input file can ');
fprintf(1,'be created.\n');
OK = FALSE;
end
end
if FLAG == 3
fprintf(1,'Input the function F(x) in terms of x.\n');
fprintf(1,'For example: sin(x)\n');
s = input(' ');
F = inline(s,'x');
fprintf(1,'Input F''(x) in terms of x.\n');
s = input(' ');
FP = inline(s,'x');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input the number of data points minus 1\n');
N = input(' ');
if N > 0
X = zeros(1,N+1);
Q = zeros(2*N+2,2*N+2);
for I1 = 0:N
fprintf(1,'Input X(%d)\n', I1);
X(I1+1) = input(' ');
Q(2*I1+1,1) = F(X(I1+1));
Q(2*I1+2,2) = FP(X(I1+1));
end
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n');
end
end
end
if OK == TRUE
% STEP 1
Z = zeros(2*N+2);
for I = 0:N
% STEP 2
Z(2*I+1) = X(I+1);
Z(2*I+2) = X(I+1);
Q(2*I+2,1) = Q(2*I+1,1);
% STEP 3
if I ~= 0
Q(2*I+1,2) = (Q(2*I+1,1)-Q(2*I,1))/(Z(2*I+1)-Z(2*I));
end
end
% STEP 4
K = 2*N+1;
for I = 2:K
for J = 2:I
Q(I+1,J+1) = (Q(I+1,J)-Q(I,J))/(Z(I+1)-Z(I-J+1));
end
end
% STEP 5
fprintf(1,'Choice of output method:\n');
fprintf(1,'1. Output to screen\n');
fprintf(1,'2. Output to text file\n');
fprintf(1,'Please enter 1 or 2\n');
FLAG = input(' ');
if FLAG == 2
fprintf(1,'Input the file name in the form - drive:\\name.ext\n');
fprintf(1,'for example: A:\\OUTPUT.DTA\n');
NAME = input(' ','s');
OUP = fopen(NAME,'wt');
else OUP = 1;
end
fprintf(OUP, 'HERMITE INTERPOLATING POLYNOMIAL\n\n');
fprintf(OUP, 'The input data follows:\n');
fprintf(OUP, ' X, F(X), F''(x)\n');
for I = 0:N
fprintf(OUP,' %12.10e %12.10e %12.10e\n',X(I+1),Q(2*I+1,1),Q(2*I+2,2));
end
fprintf(OUP, '\nThe Coefficients of the Hermite Interpolation ');
fprintf(OUP, 'Polynomial\n');
fprintf(OUP, 'in order of increasing exponent follow:\n\n');
for I = 0:K
fprintf(OUP, ' %12.10e\n', Q(I+1,I+1));
end
fprintf(1,'Do you wish to evaluate this polynomial?\n');
fprintf(1,'Enter Y or N\n');
A = input(' ','s');
if A == 'Y' | A == 'y'
fprintf(1,'Enter a point at which to evaluate\n');
XX = input(' ');
S = Q(K+1,K+1)*(XX-Z(K));
for I = 2:K
J = K-I+1;
S = (S+Q(J+1,J+1))*(XX-Z(J));
end
S = S + Q(1,1);
fprintf(OUP, 'x-value and interpolated-value\n');
fprintf(OUP, ' %12.10e %12.10e\n', XX, S);
end
if OUP ~= 1
fclose(OUP);
fprintf(1,'Output file %s created successfully\n',NAME);
end
end
% NATURAL CUBIC SPLINE ALGORITHM 3.4
%
% To construct the cubic spline interpolant S for the function f,
% defined at the numbers x(0) < x(1) < ... < x(n), satisfying
% S''(x(0)) = S''(x(n)) = 0:
%
% INPUT: n; x(0), x(1), ..., x(n); either generate A(I) = f(x(I))
% for I = 0, 1, ..., n or input A(I) for I = 0, 1, ..., n.
%
% OUTPUT: A(J), B(J), C(J), D(J) for J = 0, 1, ..., n - 1.
%
% NOTE: S(x) = A(J) + B(J)*( x - x(J) ) + C(J)*( x - x(J) )**2 +
% D(J) * ( x - x(J) )**3 for x(J) <= x < x(J + 1)
syms('OK','FLAG','N','I','X','A','AA','NAME','INP','F','M','H','XA')
syms('XL','XU','XZ','C','J','B','D','OUP','x','s');
TRUE=1;
FALSE=0;
fprintf(1,'This is the natural cubic spline interpolation.\n');
OK = FALSE;
while OK == FALSE
fprintf(1,'Choice of input method:\n');
fprintf(1,'1. Input entry by entry from keyboard\n');
fprintf(1,'2. Input data from a text file\n');
fprintf(1,'3. Generate data using a function F with nodes entered ');
fprintf(1,'from keyboard\n');
fprintf(1,'4. Generate data using a function F with nodes from ');
fprintf(1,'a text file\n');
fprintf(1,'Choose 1, 2, 3, or 4 please\n');
FLAG = input(' ');
if FLAG >= 1 & FLAG <= 4
OK = TRUE;
end
end
if FLAG == 1
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
OK = TRUE;
X = zeros(1,N+1);
A = zeros(1,N+1);
for I = 0:N
fprintf(1,'Input X(%d) and F(X(%d)) ', I, I);
fprintf(1,'on separate lines.\n');
X(I+1) = input(' ');
A(I+1) = input(' ');
end
else fprintf(1,'Number must be a positive integer\n');
end
end
end
if FLAG == 2
fprintf(1,'Has a text file been created with the data in two ');
fprintf(1,'columns ?\n');
fprintf(1,'Enter Y or N\n');
AA = input(' ','s');
if AA == 'Y' | AA == 'y'
fprintf(1,'Input the file name in the form - ');
fprintf(1,'drive:\\name.ext\n');
fprintf(1,'For example: A:\\DATA.DTA\n');
NAME = input(' ','s');
INP = fopen(NAME,'rt');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
X = zeros(1,N+1);
A = zeros(1,N+1);
for I = 0:N
X(I+1) = fscanf(INP, '%f',1);
A(I+1) = fscanf(INP, '%f',1);
end
fclose(INP);
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n');
end
end
else
fprintf(1,'Please create the input file in two column ');
fprintf(1,'form with the\n');
fprintf(1,'X values and F(X) values in the corresponding columns.\n');
fprintf(1,'The program will end so the input file can be created.\n');
OK = FALSE;
end
end
if FLAG == 3
fprintf(1,'Input the function F(x) in terms of x.\n');
fprintf(1,'For example: cos(x) \n');
s = input(' ');
F = inline(s,'x');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
X = zeros(1,N+1);
A = zeros(1,N+1);
for I = 0:N
fprintf(1,'Input X(%d)\n', I);
X(I+1) = input(' ');
A(I+1) = F(X(I+1));
end
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n');
end
end
end
if FLAG == 4
fprintf(1,'Has the text file with X-values been created?\n');
fprintf(1,'Enter Y or N\n');
AA = input(' ','s');
if AA == 'Y' | AA == 'y'
fprintf(1,'Input the file name in the form - ');
fprintf(1,'drive:\\name.ext\n');
fprintf(1,'For example: A:\\DATA.DTA\n');
NAME = input(' ','s');
INP = fopen(NAME,'rt');
fprintf(1,'Input the function F(x) in terms of x.\n');
fprintf(1,'For example: cos(x) \n');
s = input(' ');
F = inline(s,'x');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
OK = TRUE;
X = zeros(1,N+1);
A = zeros(1,N+1);
for I = 0:N
X(I+1) = fscanf(INP, '%f',1);
A(I+1) = F(X(I+1));
end
fclose(INP);
else fprintf(1,'Number must be a positive integer\n');
end
end
else
fprintf(1,'The program will end so the input file can be created\n');
OK = FALSE;
end
end
if OK == TRUE
M = N - 1;
% STEP 1
H = zeros(1,M+1);
for I = 0:M
H(I+1) = X(I+2) - X(I+1);
end
% STEP 2
% Use XA in place of ALPHA
XA = zeros(1,M+1);
for I = 1:M
XA(I+1) =
3.0*(A(I+2)*H(I)-A(I+1)*(X(I+2)-X(I))+A(I)*H(I+1))/(H(I+1)*H(I));
end
% STEP 3
% STEPs 3, 4, 5 and part of 6 solve the tridiagonal system using
% Crout reduction.
% use XL, XU, XZ in place of L, MU, Z resp.
XL = zeros(1,N+1);
XU = zeros(1,N+1);
XZ = zeros(1,N+1);
XL(1) = 1;
XU(1) = 0;
XZ(1) = 0;
% STEP 4
for I = 1:M
XL(I+1) = 2*(X(I+2)-X(I))-H(I)*XU(I);
XU(I+1) = H(I+1)/XL(I+1);
XZ(I+1) = (XA(I+1)-H(I)*XZ(I))/XL(I+1);
end
% STEP 5
XL(N+1) = 1;
XZ(N+1) = 0;
B = zeros(1,N+1);
C = zeros(1,N+1);
D = zeros(1,N+1);
C(N+1) = XZ(N+1);
% STEP 6
for I = 0:M
J = M-I;
C(J+1) = XZ(J+1)-XU(J+1)*C(J+2);
B(J+1) = (A(J+2)-A(J+1))/H(J+1) - H(J+1) * (C(J+2) + 2.0 * C(J+1)) /
3.0;
D(J+1) = (C(J+2) - C(J+1)) / (3.0 * H(J+1));
end
% STEP 7
fprintf(1,'Select output destination\n');
fprintf(1,'1. Screen\n');
fprintf(1,'2. Text file\n');
fprintf(1,'Enter 1 or 2\n');
FLAG = input(' ');
if FLAG == 2
fprintf(1,'Input the file name in the form - drive:\\name.ext\n');
fprintf(1,'For example: A:\\OUTPUT.DTA\n');
NAME = input(' ','s');
OUP = fopen(NAME,'wt');
else
OUP = 1;
end
fprintf(OUP, 'NATURAL CUBIC SPLINE INTERPOLATION\n\n');
fprintf(OUP, 'The numbers X(0), ..., X(N) are:\n');
for I = 0:N
fprintf(OUP, ' %11.8f', X(I+1));
end
fprintf(OUP, '\n\nThe coefficients of the spline on the subintervals ');
fprintf(OUP, 'are:\n');
fprintf(OUP, 'for I = 0, ..., N-1\n');
fprintf(OUP, ' A(I) B(I) C(I) D(I)\n');
for I = 0:M
fprintf(OUP,'%13.8f %13.8f %13.8f %13.8f
\n',A(I+1),B(I+1),C(I+1),D(I+1));
end
if OUP ~= 1
fclose(OUP);
fprintf(1,'Output file %s created successfully \n',NAME);
end
end
% CLAMPED CUBIC SPLINE ALGORITHM 3.5
%
% To construct the cubic spline interpolant S for the function f,
% defined at the numbers x(0) < x(1) < ... < x(n), satisfying
% S'(x(0)) = f'(x(0)) and S'(x(n)) = f'(x(n)):
%
% INPUT: n; x(0), x(1), ..., x(n); either generate A(I) = f(x(I))
% for i = 0, 1, ..., n or input A(I) for I = 0, 1, ..., n;
% FP0 = f'(x(0)); FPN = f'(x(n)).
%
% OUTPUT: A(J), B(J), C(J), D(J) for J = 0, 1, ..., n - 1.
%
% NOTE: S(x) = A(J) + B(J) * ( x - x(J) ) + C(J) * ( x - x(J) )**2 +
% D(J) * ( x - x(J) )**3 for x(J) <= x < x(J + 1)
syms('OK', 'FLAG', 'N', 'I', 'X', 'A', 'AA', 'NAME', 'INP');
syms('FP0', 'FPN', 'M', 'H', 'XA', 'XL', 'XU', 'XZ', 'C', 'J');
syms('B', 'D', 'OUP','s','x');
TRUE = 1;
FALSE = 0;
fprintf(1,'This is Clamped Cubic Spline Interpolation.\n');
OK = FALSE;
while OK == FALSE
fprintf(1,'Choice of input method:\n');
fprintf(1,'1. Input entry by entry from keyboard\n');
fprintf(1,'2. Input data from a text file\n');
fprintf(1,'3. Generate data using a function F with nodes entered ');
fprintf(1,'from keyboard\n');
fprintf(1,'4. Generate data using a function F with nodes from ');
fprintf(1,'a text file\n');
fprintf(1,'Choose 1, 2, 3, or 4 please\n');
FLAG = input(' ');
if FLAG >= 1 & FLAG <= 4
OK = TRUE;
end
end
if FLAG == 1
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
OK = TRUE;
X = zeros(1,N+1);
A = zeros(1,N+1);
for I = 0:N
fprintf(1,'Input X(%d) and F(X(%d)) ', I, I);
fprintf(1,'on separate lines\n');
X(I+1) = input(' ');
A(I+1) = input(' ');
end
else fprintf(1,'Number must be a positive integer\n');
end
end
end
if FLAG == 2
fprintf(1,'Has a text file been created with the data in two
columns ?\n');
fprintf(1,'Enter Y or N\n');
AA = input(' ','s');
if AA == 'Y' | AA == 'y'
fprintf(1,'Input the file name in the form - ');
fprintf(1,'drive:\\name.ext\n');
fprintf(1,'For example: A:\\DATA.DTA\n');
NAME = input(' ','s');
INP = fopen(NAME,'rt');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
X = zeros(1,N+1);
A = zeros(1,N+1);
for I = 0:N
X(I+1) = fscanf(INP, '%f',1);
A(I+1) = fscanf(INP, '%f',1);
end
fclose(INP);
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n');
end
end
else
fprintf(1,'Please create the input file in two column ');
fprintf(1,'form with the\n');
fprintf(1,'X values and F(X) values in the corresponding columns.\n');
fprintf(1,'The program will end so the input file can be created.\n');
OK = FALSE;
end
end
if FLAG == 3
fprintf(1,'Input the function F(x) in terms of x\n');
fprintf(1,'For example: cos(x)\n');
s = input(' ');
F = inline(s,'x');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
X = zeros(1,N+1);
A = zeros(1,N+1);
for I = 0:N
fprintf(1,'Input X(%d)\n', I);
X(I+1) = input(' ');
A(I+1) = F(X(I+1));
end
OK = TRUE;
else
fprintf(1,'Number must be a positive integer\n');
end
end
end
if FLAG == 4
fprintf(1,'Has the text file with X-values been created?\n');
fprintf(1,'Enter Y or N\n');
AA = input(' ','s');
if AA == 'Y' | AA == 'y'
fprintf(1,'Input the file name in the form - ');
fprintf(1,'drive:\\name.ext\n');
fprintf(1,'For example: A:\\DATA.DTA\n');
NAME = input(' ','s');
INP = fopen(NAME,'rt');
fprintf(1,'Input the function F(x) in terms of x\n');
fprintf(1,'For example: cos(x)\n');
s = input(' ');
F = inline(s,'x');
OK = FALSE;
while OK == FALSE
fprintf(1,'Input n\n');
N = input(' ');
if N > 0
OK = TRUE;
X = zeros(1,N+1);
A = zeros(1,N+1);
for I = 0:N
X(I+1) = fscanf(INP, '%f',1);
A(I+1) = F(X(I+1));
end
fclose(INP);
else fprintf(1,'Number must be a positive integer\n');
end
end
else
fprintf(1,'The program will end so the input file can be created\n');
OK = FALSE;
end
end
if OK == TRUE
fprintf(1,'Enter F''(X(0)) and F''(X(N)) on separate lines\n');
FP0 = input(' ');
FPN = input(' ');
end
if OK == TRUE
% STEP 1
M = N - 1;
H = zeros(1,M+1);
for I = 0:M
H(I+1) = X(I+2) - X(I+1);
end
% STEP 2
% use XA instead of alpha
XA = zeros(1,N+1);
XA(1) = 3.0 * (A(2) - A(1)) / H(1) - 3.0 * FP0;
XA(N+1) = 3.0 * FPN - 3.0 * (A(N+1) - A(N)) / H(N);
% STEP 3
for I = 1:M
XA(I+1) =
3.0*(A(I+2)*H(I)-A(I+1)*(X(I+2)-X(I))+A(I)*H(I+1))/(H(I+1)*H(I));
end
% STEP 4
% STEPS 4, 5, 6 and part of 7 solve the tridiagonal system using
Algorithm 6.7
% use XL, XU, XZ in place of L, MU, Z resp.
XL = zeros(1,N+1);
XU = zeros(1,N+1);
XZ = zeros(1,N+1);
XL(1) = 2.0 * H(1);
XU(1) = 0.5;
XZ(1) = XA(1) / XL(1);
% STEP 5
for I = 1:M
XL(I+1) = 2.0 * (X(I+2) - X(I)) - H(I) * XU(I);
XU(I+1) = H(I+1) / XL(I+1);
XZ(I+1) = (XA(I+1) - H(I) * XZ(I)) / XL(I+1);
end
% STEP 6
XL(N+1) = H(N) * (2.0 - XU(N));
XZ(N+1) = (XA(N+1) - H(N) * XZ(N)) / XL(N+1);
C = zeros(1,N+1);
B = zeros(1,N+1);
D = zeros(1,N+1);
C(N+1) = XZ(N+1);
% STEP 7
for I = 1:N
J = N - I;
C(J+1) = XZ(J+1) - XU(J+1) * C(J+2);
B(J+1) = (A(J+2)-A(J+1))/H(J+1)-H(J+1)*(C(J+2)+2.0*C(J+1))/3.0;
D(J+1) = (C(J+2) - C(J+1)) / (3.0 * H(J+1));
end
% STEP 8
fprintf(1,'Select output destination\n');
fprintf(1,'1. Screen\n');
fprintf(1,'2. Text file\n');
fprintf(1,'Enter 1 or 2\n');
FLAG = input(' ');
if FLAG == 2
fprintf(1,'Input the file name in the form - drive:\\name.ext\n');
fprintf(1,'For example: A:\\OUTPUT.DTA\n');
NAME = input(' ','s');
OUP = fopen(NAME,'wt');
else
OUP = 1;
end
fprintf(OUP, 'CLAMPED CUBIC SPLINE INTERPOLATION\n\n');
fprintf(OUP, 'The numbers X(0), ..., X(N) are:\n');
for I = 0:N
fprintf(OUP, ' %11.8f', X(I+1));
end
fprintf(OUP, '\n\nThe coefficients of the spline on the subintervals ');
fprintf(OUP, 'are:\n');
fprintf(OUP, 'for I = 0, ..., N-1\n');
fprintf(OUP, ' A(I) B(I) C(I) D(I)\n');
for I = 0:M
fprintf(OUP,'%13.8f %13.8f %13.8f
%13.8f\n',A(I+1),B(I+1),C(I+1),D(I+1));
end
if OUP ~= 1
fclose(OUP);
fprintf(1,'Output file %s created successfully\n',NAME);
end
end
Spline cubico gráfica
%interpcub
%programa para efectuar interpolacion cubica por segmentos para datos
% no igualmente espaciados
clear all
%Entrada de informacion
n=input('numero de datos=');
x=input('vector de abscisas=');
y=input('vector de ordenadas=');
W1=input('segunda derivada al principio=');
WN=input('segunda derivada al final=');
%Calculo de los tamaños de paso
for i=1:n-1
h(i)=x(i+1)-x(i);
end
%calculo de las segundas diferencias divididas
for k=1:n-2
b(k)=6*((y(k+2)-y(k+1))/h(k+1)-(y(k+1)-y(k))/h(k));
end
%Sistema de ecuaciones para hallar las segundas derivadas, desde w(2)
hasta w(n-1)
% M*W=B
M(1,1)=2*(h(1)+h(2));
M(1,2)=h(2);
B(1)=b(1)-h(1)*W1;
M(n-2,n-3)=h(n-2);
M(n-2,n-2)=2*(h(n-2)+h(n-1));
B(n-2)=b(n-2)-h(n-1)*WN;
for k=2:n-3
M(k,k-1)=h(k);
M(k,k)=2*(h(k)+h(k+1));
M(k,k+1)=h(k+1);
B(k)=b(k);
end
aux=inv(M)*B';
W(1)=W1;
W(n)=WN;
for j=2:n-1;
W(j)=aux(j-1);
end
%Calculo de las primeras derivadas
for k=1:n-1
U(k)=(y(k+1)-y(k))/h(k)-(h(k)/6)*(2*W(k)+W(k+1));
end
% Calculo de los coeficientes de los polinomios
for k=1:n-1
c(k,1)=(W(k+1)-W(k))/(6*h(k));
c(k,2)=W(k)/2-3*x(k)*c(k,1);
c(k,3)=U(k)-x(k)*W(k)+3*c(k,1)*x(k)^2;
c(k,4)=-c(k,1)*x(k)^3+(x(k)^2)*W(k)/2-x(k)*U(k)+y(k);
end
%Representacion grafica de los polinomios
for i=1:n-1
t(i,:)=x(i):0.01*h(i):x(i+1);
p(i,:)=c(i,1).*(t(i,:).^3)+c(i,2).*(t(i,:).^2)+c(i,3).*t(i,:)+c(i,4);
plot(t(i,:),p(i,:),'k')
grid on
hold on
end
plot (x,y,'o')
for j=1:n-1
for i=1:4
s(i)=c(j,i);
end
disp(strcat('p',num2str(j),'=',poly2str(s,'t')));
end
