Programa 5.3 (Cercha cubica sujeta - pagina 316)

%Programa 5.3 (Cercha cubica sujeta - pagina 316)
% Construccion y determinacion de la cercha cubica sujeta S(x) que interpola los N + 1
% datos {(Xk, Yk)}, con k=1,2,3...N
%
%
%Sintaxis
% csfit(X,Y,dx0,dxn)
%
%Entrada
% - X es el vector fila de abscisas 1 x n
% - Y es el vector fila de ordenadas 1 x n
% - dx0 = S'(x0) es la derivada en el primer extremo
% - dxn = S'(xn) es la derivada en el segundo extremo
%Resultados
% - S las filas de S son los coeficientes, en orden decreciente
% en cada cubica de la cercha

N=length(X)-1;
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);
C=H(2:N);
U=6*diff(D);

%Clamped spline endpoint constraints

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));

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

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

%Clamped spline endpoint constraints

M(1)=3*(D(1)-dx0)/H(1)-M(2)/2;
M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2;

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

Programa 5.2 (Polinoio optimo en minimos cuadrados - pagina 296).

%Programa 5.2 (Polinoio optimo en minimos cuadrados - pagina 296).
% Construccion del polinomio de grado M dado por
% Pm(x)=C1 + C2·X + C3·X^2 + ··· + Cn·X^M
% que mejor se ajusta en el sentido de los minimos cuadrados a los N datos (X1,Y1),...,(Xn,Yn).
%
%
%Sintaxis
% lspoly(X,Y,M)
%
%Entrada
% - X es el vector fila de abscisas 1 x n
% - Y es el vector fila de ordenadas 1 x n
% - M es el grado del polinomio optimo
%Resultados
% - C es el vector de coeficientes del polinomio
% en potencias decrecientes de x

n=length(X);
B=zeros(1:M+1);
F=zeros(n,M+1);

%Fill the columns of F with the powers of X

for k=1:M+1
F(:,k)=X'.^(k-1);
end

%Solve the linear system from (25)

A=F'*F;
B=F'*Y';
C=A\B;
C=flipud(C);

Programa 5.1 (Recta de regresion - pagina 280).

%Programa 5.1 (Recta de regresion - pagina 280).
% Construccion de la recta de regresion y = Ax + B que mejor ajusta en el sentido de los
% minimos cuadrados a los N datos (X1,Y1),...,(Xn,Yn).
%
% Se resuelve el sig sistema para obtener los coef. A y B de la recta que mejor ajusta.
% suma(Xk^2) · A + suma(Xk) · B = suma(Xk·Yk)
% suma(Xk) · A + N · B = suma(Yk)
%
%Sintaxis
% lsline(X,Y)
%
%Entrada
% - X es el vector fila de abscisas 1 x n
% - Y es el vector fila de ordenadas 1 x n
%Resultados
% - A es el coeficiente de x en Ax + B
% - B es el termino independiente en Ax + B

xmean=mean(X);
ymean=mean(Y);
sumx2=(X-xmean)*(X-xmean)';
sumxy=(Y-ymean)*(X-xmean)';
A=sumxy/sumx2;
B=ymean-A*xmean;

Programa 4.2 (Polinomio interpolador de Newton - pagina 247).

%Programa 4.2 (Polinomio interpolador de Newton - pagina 247).
% Construccion del polinomio interpolador de Newton de grado menor o igual que N que
% pasa por los puntos (Xk,Yk)=(Xk,f(Xk)) para k=0,1,...N.
% P(x)=d<0,0> + d<1,1>·(X-X<0>) + d<2,2>·(X-X<0>)·(X-X<1>)+···
% d·(X-X<0>)·(X-X<1>)···(X-X)
% siendo
% d=y y d=(d - d) / (d - d)
%
%Sintaxis
% [P,D]=newpoly(X,Y)
%
%Entrada
% - X es un vector que contiene la lista de las abscisas.
% - Y es un vector que contiene la lista de las ardenadas.
%Resultados
% - P es la matriz que contiene los coeficientes del polinomio interpolador
% - de Newton, escrito en forma habitual, en potencias decrecientes de x.
% - D es la tabla de diferencias divididas, para construir el polinomio se toman los elementos de la
% diagonal utilizando la siguiente formula (d·(X-X<0>)·(X-X<1>)···(X-X)

n=length(X);
D=zeros(n,n);
D(:,1)=Y';

%Usamos la formula d=(d - d) / (d - d)
%para hallar la tabla de diferencias divididas

for j=2:n
for k=j:n
D(k,j)=(D(k,j-1)-D(k-1,j-1))/(X(k)-X(k-j+1));
end
end

%Calculo del vector que contiene los coeficientes del polinomio interpolador
%de Newton, escrito en forma habitual, en potencias decrecientes de x.

P=D(n,n);

for k=(n-1):-1:1
P=conv(P,poly(X(k)));
m=length(P);
P(m)=P(m)+D(k,k);
end

Programa 4.1 (Polinomio interpolador de Lagrange - pagina 235).

%Programa 4.1 (Polinomio interpolador de Lagrange - pagina 235).
% Construccion del polinomio interpolador de Lagrange
% P(x) = sumatoria(Y·L(x)) , con k=0,1,...N
% que pasa por los N+1 puntos (X, Y) para k=0,1,...N.
% y L=Multiplicatoria(X - X) / Multiplicatoria(X - X) siempre que X~=X
% j=0,1,2,...N
%
%Sintaxis
% [P,L]=lagran(X,Y)
%
%Entrada
% - X es un vector que contiene la lista de las abscisas.
% - Y es un vector que contiene la lista de las ardenadas.
%Resultados
% - P es la matriz que contiene los coeficientes del
% - polinomio interpolador de Lagrange en orden decreciente.
% - L es la matriz que contiene los coeficientes de los
% - polinomios coeficientes de Lagrange en odren decreciente.

w=length(X);
n=w-1;
L=zeros(w,w);

%Formacion de los polinomios coeficientes de Lagrange.

for k=1:n+1
V=1;
for j=1:n+1
if k~=j
V=conv(V,poly(X(j)))/(X(k)-X(j));
end
end
L(k,:)=V;
end

%Calculo de los coeficientes del polinomio
%interpolador de Lagrange

P=Y*L

Programa 3.5 (Metodo iterativo de Gauss-Seidel - pagina 180).

%Programa 3.5 (Metodo iterativo de Gauss-Seidel - pagina 180).
% Resolucion e un sistema lineal AX=B mediante la generacion de una sucesion {P}
% que converge a la solucion, a partir de un punto inicial P<0>. Una condicion
% suficiente para que el metodo sea aplicable es que A sea diagonal estrictamente dominante.
%
%Sintaxis
% X=gseid(A,B,P,delta, max1)
%
%Entrada
% - A es una matriz invertible de orden de orden n x n (cuadrada)
% - B es una matriz de orden n x 1 (columna)
% - P es una matriz de orden N x 1: el punto inicial
% - delta es la tolerancia de P
% - max1 es el numero maximo de iteraciones
%Resultados
% - X es la matriz de orden N x 1:
% la aproximacion a la solucion de AX=B
% generada por el metodo iterativo de Gauss-Seidel

N = length(B);

for k=1:max1
for j=1:N
if j==1
X(1)=(B(1)-A(1,2:N)*P(2:N))/A(1,1);
elseif j==N
X(N)=(B(N)-A(N,1:N-1)*(X(1:N-1))')/A(N,N);
else
%X contains the kth approximations and P the (k-1)st
X(j)=(B(j)-A(j,1:j-1)*X(1:j-1)-A(j,j+1:N)*P(j+1:N))/A(j,j);
end
end
err=abs(norm(X'-P));
relerr=err/(norm(X)+eps);
P=X';
if (err < delta) (relerr < delta)
break
end
end

X=X';

Programa 3.4 (Metodo iterativo de Jacobi - pagina 179).

%Programa 3.4 (Metodo iterativo de Jacobi - pagina 179).
% Resolucion e un sistema lineal AX=B mediante la generacion de una sucesion {P}
% que converge a la solucion, a partir de un punto inicial P<0>. Una condicion
% suficiente para que el metodo sea aplicable es que A sea diagonal estrictamente dominante.
%
%Sintaxis
% X=jacobi(A,B,P,delta, max1)
%
%Entrada
% - A es una matriz invertible de orden n x n (cuadrada)
% - B es una matriz de orden n x 1 (columna)
% - P es una matriz de orden N x 1: el punto inicial
% - delta es la tolerancia de P
% - max1 es el numero maximo de iteraciones
%Resultados
% - X es la matriz de orden N x 1:
% la aproximacion a la solucion de AX=B
% generada por el metodo iterativo de Jacobi

N = length(B);

for k=1:max1
for j=1:N
X(j)=(B(j)-A(j,[1:j-1,j+1:N])*P([1:j-1,j+1:N]))/A(j,j);
end
err=abs(norm(X'-P));
relerr=err/(norm(X)+eps);
P=X';
if (err < delta) (relerr < delta)
break
end
end

X=X';