Programa 7.3 (Regla recursiva del trapecio - pagina 409).

%Programa 7.3 (Regla recursiva del trapecio - pagina 409).
% Construccion de la aproximacion a la integral
%
% b 2^j
% I f(x)dx = h/2 Sum (f(x) + f(x))
% a k=1
%
% usando recursivamente la regla del trapecio conforme se incrementa el numero
% de subintervalos de [a,b]. En la iteracion J-esima se toman los valores de f(x)
% en 2^j + 1 nodos equiespaciados
%
%Sintaxis
% rctrap(f,a,b,n)
%
%Datos
% - f es el integrando, dado como una cadena de caracteres 'f'
% - a y b son los extremos inferior y superior del
% intervalo de integracion
% - n es el numero de veces que se hce la recursion
%Resultados
% - T es la lista de las aproximaciones obtenidas con la regla
% recursiva del trapecio

M=1;
h=b-a;
T=zeros(1,n+1);
x=a;
fa=eval(f);
x=b;
fb=eval(f);
T(1)=h*(fa+fb)/2;

for j=1:n
M=2*M;
h=h/2;
s=0;
for k=1:M/2
x=a+h*(2*k-1);
s=s+eval(f);
end
T(j+1)=T(j)/2+h*s;
end

Programa 7.2 (Regla Compuest de Simpson - pagina 393).

%Programa 7.2 (Regla Compuest de Simpson - pagina 393).
% Construccion de la aproximacion a la integral
%
% b M-1 M
% I f(x)dx = h/3 (f(a)+f(b)) + 2h/3 Sum f(x<2k>)+ 4h/3 Sum f(x<2k-1>)
% a k=1 k=1
%
% evaluando f(x) en los 2M + 1 nodos equiespaciados x=a+kh,
% para k=0,1,2...,2M. Notese que X<0>=a y que x<2m>=b.
%
%Sintaxis
% simprl(f,a,b,M)
%
%Datos
% - f es el integrando, dado como una cadena de caracteres 'f'
% - a y b son los extremos inferior y superior del
% intervalo de integracion
% - M es el numero de suintervalos
%Resultados
% - s es la aproximacion obtenida con la regla compuesta del Simpson

h=(b-a)/(2*M);
s1=0;
s2=0;

for k=1:M
x=a+h*(2*k-1);
s1=s1+eval(f);
end
for k=1:(M-1)
x=a+h*2*k;
s2=s2+eval(f);
end
x=a;
fa=eval(f);
x=b;
fb=eval(f);
s=h*(fa+fb+4*s1+2*s2)/3;

Programa 7.1 (Regla Compuest del trapecio - pagina 393).

%Programa 7.1 (Regla Compuest del trapecio - pagina 393).
% Construccion de la aproximacion a la integral
%
% b M-1
% I f(x)dx = h/2 (f(a)+f(b)) + h Sum f(x)
% a k=1
%
% evaluando f(x) en los M + 1 nodos equiespaciados x=a+kh,
% para k=0,1,2...,M. Notese que X<0>=a y que x=b.
%
%Sintaxis
% traprl(f,a,b,M)
%
%Datos
% - f es el integrando, dado como una cadena de caracteres 'f'
% - a y b son los extremos inferior y superior del
% intervalo de integracion
% - M es el numero de suintervalos
%Resultados
% - s es la aproximacion abtenida con la regla compuesta del trapecio


h=(b-a)/M;
s=0;

for k=1:(M-1)
x=a+h*k;
s=s+eval(f);
end
x=a;
fa=eval(f);
x=b;
fb=eval(f);
s=h*(fa+fb)/2+h*s;

Programa 6.3 (Derivacion basada en N + 1 nodos - pagina 365).

%Programa 6.3 (Derivacion basada en N + 1 nodos - pagina 365).
% Construccion del polinomio interpolador de Newton de grado N
% P(x)=d<0> + d<1>·(X-X<0>) + d<2>·(X-X<0>)·(X-X<1>)+···
% d·(X-X<0>)·(X-X<1>)···(X-X)
% para aproximar numericamente f'(x) usando f'(x<0>)=P'(x<0>) como respuesta final.
% El metodo debe usarse solo con x<0>. Los nodos pueden reordenarse como
% {X,X<0>,...,X,X,...,X} si se desa calcular f'(X)=P'(X)
%
%Sintaxis
% [A,df]=diffnew(X,Y)
%
%Entrada
% - X es un vector de 1 x n que contiene las absisas
% - Y es un vector de 1 x n que contiene las ordenadas
%Resultados:
% - A es el resultado 1 x n que contiene los coeficientes del
% polinomio de Newton de grado N
% - df es la derivada aproximada

A=Y;
N=length(X);

for j=2:N
for k=N:-1:j
A(k)=(A(k)-A(k-1))/(X(k)-X(k-j+1));
end
end

x0=X(1);
df=A(2);
prod=1;
n1=length(A)-1;

for k=2:n1
prod=prod*(x0-X(k));
df=df+prod*A(k+1);
end

Programa 6.2 (Derivacion numerica mediante extrapolacion - pagina 349).

%Programa 6.2 (Derivacion numerica mediante extrapolacion - pagina 349).
% Construccion de la tabla D(j,k) (con k<=j) de aproximaciones numericas a f'(x) en la que se
% utiliza f'(x)=D(n,n) como respuesta final. Las aproximaciones D(j,k) se almacenan en una
% matriz triangular inferior cuya primera columna viene dada por
%
% D(j,k) = ( ( f(x+c) - f(x-c) ) / (2·c) ), con c=2^(-j)·h
% y cuyos elementos en la fila j-esima, para j>=2, son
% D(j,k) = D(j,k-1) + ( D(j,k-1) - D(j-1,k-1) ) / (4^k - 1), donde 2<=k<=j.
%
%Sintaxis
% [D,err,relerr,n]=diffext(f,x,delta,toler)
%
%Entrada
% - f es la funcion, introducida como una cadena de caracteres 'f'
% - x es el punto en el que se deriva
% - delta es la toleracia para el error
% - toler es la tolerancia para el error relativo
%Resultados:
% - D es la matriz de las aproximaciones de la derivada
% - err es la cota de error
% - relerr es la cota de error relativo
% - n es la coordenada de la "mejor aproximacion"

err=1;
relerr=1;
h=1;
j=1;
D(1,1)=(feval(f,x+h)-feval(f,x-h))/(2*h);

while relerr > toler & err > delta &j <12
h=h/2;
D(j+1,1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
for k=1:j
D(j+1,k+1)=D(j+1,k)+(D(j+1,k)-D(j,k))/((4^k)-1);
end
err=abs(D(j+1,j+1)-D(j,j));
relerr=2*err/(abs(D(j+1,j+1))+abs(D(j,j))+eps);
j=j+1;
end

[n,n]=size(D);

Programa 6.1 (Derivacion numerica mediante limites - pagina 348).

%Programa 6.1 (Derivacion numerica mediante limites - pagina 348).
% Construccion de las aproximaciones numericas a f'(x) mediante la generacion de
% una sucesion
% f'(x)=D=( ( f(x+c) - f(x-c) / (2·c) ), con c=10^(-k)·h y k=0,1,2,...,n
%que continua hasta que
% ( D - D >= D - D ) o bien ( toler > D - D )
% que es el criterio que se trata de encontrar para la mejor aproximacion D=f'(x).
%
%Sintaxis
% [L,n]=difflim(f,x,toler)
%
%Entrada
% - f es la funcion, introducida como una cadena de caracteres 'f'
% - x es el punto en el que se deriva
% - toler es la tolerancia para el error
%Resultados:
% - L=[ h f=f(1+h) f-e ], donde:
% H es el vector de los incrementos (
% D es el vector de las aproximaciones a la derivada
% E es el vector de las cotas del error
% Al construir la tabla de diferencias incrementales se adiciona una columna mas D = (f-e) / h
% - n es la coordenada de la "mejor aproximacion"

max1=15;
h=1;
H(1)=h;
D(1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
E(1)=0;
R(1)=0;

for n=1:2
h=h/10;
H(n+1)=h;
D(n+1)=(feval(f,x+h)-feval(f,x-h))/(2*h);
E(n+1)=abs(D(n+1)-D(n));
R(n+1)=2*E(n+1)*(abs(D(n+1))+abs(D(n))+eps);
end

n=2;

while((E(n) > E(n+1)) & (R(n) > toler)) & n < max1
h=h/10;
H(n+2)=h;
D(n+2)=(feval(f,x+h)-feval(f,x-h))/(2*h);
E(n+2)=abs(D(n+2)-D(n+1));
R(n+2)=2*E(n+2)*(abs(D(n+2))+abs(D(n+1))+eps);
n=n+1;
end

n=length(D)-1;
L=[H' D' E'];

Programa 5.4 (Polinomios Trigonometricos - pagina 331)

%Programa 5.4 (Polinomios Trigonometricos - pagina 331)
% Construccion del polinomio la cercha cubica sujeta S(x) que interpola los N + 1
% datos {(Xk, Yk)}, con k=1,2trigonometrico de grado M de la forma
% P(x) = A0/2 + sumatoria( (Aj · cos(j·X)) + (Bj · sen(j·X)) )
% correspondiente a N nodos equiespaciados Xk=-pi + (2·pi·k) / N, para k=1,2,...N.
% Esta construccion puede efectuarse si 2·M + 1 es menor o ingual a N.
%
%Sintaxis
% tpcoeff(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 trigonometrico
%Resultados
% - A es el vector de coeficientes de los cos(j·X)
% - B es el vector de coeficientes de los sen(j·X)

N=length(X)-1;
max1=fix((N-1)/2);

if M>max1
M=max1;
end

A=zeros(1,M+1);
B=zeros(1,M+1);
Yends=(Y(1)+Y(N+1))/2;
Y(1)=Yends;
Y(N+1)=Yends;
A(1)=sum(Y);

for j=1:M
A(j+1)=cos(j*X)*Y';
B(j+1)=sin(j*X)*Y';
end

A=2*A/N;
B=2*B/N;
A(1)=A(1)/2;

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

Programa 3.3 (PA=LU: factorizacion con pivoteo - pagina 166).

%Programa 3.3 (PA=LU: factorizacion con pivoteo - pagina 166).
% Calculo de la solucion del sistema lineal AX=B cuando la matriz es invertible.
%
%Sintaxis
% X = lufact(A,B)
%
%Entrada
% - A es una matriz de orden n x n (cuadrada)
% - B es una matriz de orden n x 1 (columna)
%Resultados
% - X es la matriz de orden N x 1 que contiene la solucin de AX=B.

%Inicializamos X, Y, la matriz de almacenamiento temporal C y
%la matriz fila R donde se registran los intercambios de filas
[N,N]=size(A);
X=zeros(N,1);
Y=zeros(N,1);
C=zeros(1,N);
R=1:N;

for p=1:N-1

%Determinacion de la fila pivote para la columna q-esima
[max1,j]=max(abs(A(p:N,p)));

%Intercambio de las filas q-esima y j-esima
C=A(p,:);
A(p,:)=A(j+p-1,:);
A(j+p-1,:)=C;
d=R(p);
R(p)=R(j+p-1);
R(j+p-1)=d;

if A(p,p)==0
'A es singular, no hay solucion o no es unica'
break
end

%Calculo del multiplicador,
%que se guarda en la parte subdiagonal de A
for k=p+1:N
mult=A(k,p)/A(p,p);
A(k,p) = mult;
A(k,p+1:N)=A(k,p+1:N)-mult*A(p,p+1:N);
end
end

%Resolucion para hallar Y
Y(1) = B(R(1));
for k=2:N
Y(k)= B(R(k))-A(k,1:k-1)*Y(1:k-1);
end

%Resolucion para hallar X
X(N)=Y(N)/A(N,N);
for k=N-1:-1:1
X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);
end

Programa 3.2 (Triangularizacion superior seguida de sustitucion regresiva - pagina 150).

%Programa 3.2 (Triangularizacion superior seguida de sustitucion regresiva - pagina 150).
% Calculo de la solucion del sistema lineal AX=B mediante la reduccion a la forma
% triangular superior de la matriz ampliada [AB] seguida de la sustitucion regresiva.
%
%Sintaxis
% X = uptrbk(A,B)
%
%Entrada
% - A es una matriz invertible de orden n x n (cuadrada)
% - B es una matriz de orden n x 1 (columna)
%Resultados
% - X es la matriz de orden N x 1 que contiene la solucin de AX=B.

%Inicializamos X con una matriz C que sirve de almacen temporal
[N N]=size(A);
X=zeros(N,1);
C=zeros(1,N+1);

%Calculo de la matriz ampliada Aug=[AB]
Aug=[A B];

for p=1:N-1
%Pivoteo parcial en la columna q-esima
[Y,j]=max(abs(Aug(p:N,p)));
%Intercambio las filas q-esima y (j+q-1)-esima
C=Aug(p,:);
Aug(p,:)=Aug(j+p-1,:);
Aug(j+p-1,:)=C;

if Aug(p,p)==0
'A es singular. No hay solucion o no es unica'
break
end
%Proceso de eliminacion en la columna q-esima
for k=p+1:N
m=Aug(k,p)/Aug(p,p);
Aug(k,p:N+1)=Aug(k,p:N+1)-m*Aug(p,p:N+1);
end
end

%Sustitucion regresiva en [UY] usando el programa 3.1
X=backsub(Aug(1:N,1:N),Aug(1:N,N+1));

Programa 3.1 (Sustitucion regresiva - pagina 135).

%Programa 3.1 (Sustitucion regresiva - pagina 135).
% Resolucion de un sistema triangular superior AX=B por el metodo de
% sustitucion regresiva. El metodo funciona solo si todos los elementos
% diagonales son distintos de cero. Primero se calcula X=b / a
% y luego se usa la regla
% X=( b - sumatoria(a·X) )/a ,para j=k+1,...,N y k=N-1,N-2,...,1.
%
%Sintaxis
% X=backsub(A,B)
%
%Entrada
% - A es una matriz triangular superior invertible de orden n x n (cuadrada)
% - B es una matriz de orden n x 1 (columna)
%Resultados
% - X es la solucin del sistema lineal

%Calculo de la dimencion de B e inicializacion de X
n=length(B);
X=zeros(n,1);
X(n)=B(n)/A(n,n);

for k=n-1:-1:1
X(k)=(B(k)-A(k,k+1:n)*X(k+1:n))/A(k,k);
end

Programa 2.6 (Metodo de la secante - pagina 93).

%Programa 2.6 (Metodo de la secante - pagina 93).
% Aproximacion a una raiz de f(x)=0 a partir de los valores iniciales
% p<0> y p<1> mediante la iteracion
% p = p - ( f(p)(p - p) ) / ( f(p) - f(p) ) para k=1,2,...
%
%
%Sintaxis
% secant(f,p0,p1,delta,epsilon,max1)
%
%Datos
% - f es la funcion, introducida como una cadena de caracteres 'f'
% - p0 y p1 son las aproximaciones iniciales a un cero de f
% - delta es la tolerancia para p1
% - epsilon es la tolerancia para los valores de la funcion
% - max1 es el numero maximo de iteraciones
%Resultados
% - p1 es la aproximacion al cero, optenida con el metodo de la secante
% - err es una etimacion del error de p1
% - k es el numero de iteraciones realizadas
% - y es el valor de la funcion f(p1)

for k=1:max1
x=p0;
fenp0=feval(f,x);
x=p1;
fenp1=feval(f,x);
p2=p1-fenp1*(p1-p0)/(fenp1-fenp0);
disp(['p(',num2str(k+1),')=']);disp(p2);
err=abs(p2-p1);
relerr=2*err/(abs(p2)+delta);
p0=p1;
p1=p2;
x=p1;
y=feval(f,x);
%disp(['y(',num2str(k+1),')=']);disp(y);
if (err < delta) (relerr < delta) (abs(y) < epsilon),break,end
end

Programa 2.5 (Iteracion de Newton-Raphson - pagina 92).

%Programa 2.5 (Iteracion de Newton-Raphson - pagina 92).
% Aproximacion a una raiz de f(x)=0 a partir de un valor inicia p<0>
% mediante la iteracion
% p = p - ( f(p) / f'(p) ) para k=1,2,...
%
%
%Sintaxis
% newton(f,df,p0,delta,epsilon,max1)
%
%Entrada
% - f es la funcion, introducida como una cadena de caracteres 'f'
% - df es la derivada de f, introducida como una cadena 'df'
% - p0 es la aproximacion inicial a un cero de f
% - delta es la tolerancia para p0
% - epsilon es la tolerancia para los valores de la funcion
% - max1 es el numero maximo de iteraciones
%Resultados
% - p0 es la aproximacion al cero, optenida con el metodo de Newton-Raphson
% - err es una estimacion dl error de p0
% - k es el numero de estimaciones realizadas
% - y es el valor de la funcion f(p0)


x=p0
for k=1:max1
p1=p0-feval(f,p0)/feval(df,p0);
disp(['p(',num2str(k),')=']);disp(p1);
err=abs(p1-p0);
relerr=2*err/(abs(p1)+delta);
p0=p1;
x=p0;
y=feval(f,p0)
if (err < delta) (relerr < delta) (abs(y) < epsilon),break,end
end

Programa 2.4 (Localizacion aproximada de raices - pagina 75).

%Programa 2.4 (Localizacion aproximada de raices - pagina 75).
% Estimacion aproximada de la localizacion de una raiz de la ecuacion
% f(x)=0 en el intervalo [a,b] mediante el uso de puntos de muestra equiespaciados
% (x,f(x)) y de acuerdo con los siguientes criterios:
% i) (y)(y) < 0, o
% ii) y <> - y)(y - y) < 0
% Esto es, o bien f(x) y f(x) tienen distinto signo, o bien f(x) es pequeño
% y la pendiente de la curva y=f(x) cambia de signo cerca de (x,f(x)).
%
%
%Sintaxis
% approot (f,x,epsilon)
%
%Entrada
% - f es la funcion, entrada como cadena de caracteres
% - x es el vector de abscisas
% - epsilon es la tolerancia
%Resultados
% - R es el vector de raices aproximadas


Y=eval(f);
yrange = max(Y)-min(Y);
epsilon2 = yrange*epsilon;
n=length(x);
m=0;
x(n+1)=x(n);
Y(n+1)=Y(n);

for k=2:n,
if Y(k-1)*Y(k) <= 0,
m=m+1;
R(m)=(x(k-1)+x(k))/2
end
s=(Y(k)-Y(k-1))*(Y(k+1)-Y(k));
if (abs(Y(k)) < epsilon2) & (s <= 0),
m=m+1;
R(m)=x(k)
end
end

Programa 2.3 (Metodo de la regula falsi o de la posicion falsa - pagina 66).

Programa 2.3 (Metodo de la regula falsi o de la posicion falsa - pagina 66).
% Aproximacion a una raiz de la ecuacion f(x)=0 en el intervalo [a,b].
% Puede usarse solo si f(x) es continua y f(a) y f(b) tienen distinto signo.
%
%
%Sintaxis
% regula(f,a,b,delta,epsilon,max1)
%
%Entrada
% - f es la funcion, introducida como una cadena de caracteres 'f'
% - a y b son el extremo izquierdo y el extremo derecho
% - delta es la tolerancia
% - epsilon es la tolerancia para el valor de f en el cero
% - max1 es el numero maximo de iteraciones
%Resultados
% - c es el cero
% - yc=f(c)
% - err es el error estimado de la aproximacion a c


x=a
ya=eval(f)
x=b
yb=eval(f)

if ya*yb>0
disp('Note: f(a)*f(b) >0'),
break,
end

for k=1:max1
dx=yb*(b-a)/(yb-ya);
c=b-dx;
ac=c-a;
x=c
yc=eval(f)
if yc==0,break;
elseif yb*yc > 0
b=c;
yb=yc;
else
a=c;
ya=yc;
end
dx=min(abs(dx),ac);
if abs(dx) < delta,break,end
if abs(yc) < epsilon, break,end
end

Programa 2.2 (Metodo de biseccion - pagina 65).

%Programa 2.2 (Metodo de biseccion - pagina 65). Aproximacion a una raiz de la
% ecuacion f(x)=0 en el intervalo [a,b]. Puede usarse solo si f(x) es continua
% y f(a) y f(b) tienen distinto signo.
%
%
%Sintaxis
% bisect(f,a,b,delta)
%
%Datos
% - f es la funcion, introducida como una cadena de caracteres 'f'
% - a y b son el extremo izquierdo y el extremo derecho
% - delta es la tolerancia
%Resultados
% - c es el cero
% - yc=f(c)
% - err es el error estimado de la aproximacion a c


x=a
ya=eval(f)
x=b
yb=eval(f)
if ya*yb > 0,break,end
max1=1+round((log(b-a)-log(delta))/log(2));
for k=1:max1
c=(a+b)/2;
x=c
yc=eval(f)
if yc==0
a=c;
b=c;
elseif yb*yc>0
b=c;
yb=yc;
else
a=c;
ya=yc;
end
if b-a < delta, break,end
end

c=(a+b)/2;
err=abs(b-a);
x=c
yc=eval(f)

Programa 2.1 (Iteracion de punto fijo - pagina 54).

%Programa 2.1 (Iteracion de punto fijo - pagina 54). Aproximacion a una solucion
% de la ecuacion x=g(x) mediante la iteracion p=g(p) realizada a partir
% de una aproximacin inicial p<0>
%
% Datos
% - g es la funcion de iteracion
% - p0 es el punto de partida
% - tol es la tolerancia
% - max1 es el numero maximo de iteraciones
% Resultados
% - k es el numero de iteraciones realizadas
% - p es la aproximacion al punto fijo
% - err es la diferencia entre dos terminos consecutivos
% - P el la sucecion {pn} completa

P(1)= p0;

for k=2:max1
P(k)=feval(g,P(k-1));
err=abs(P(k)-P(k-1));
relerr=err/(abs(P(k))+eps);
p=P(k);
if (err < tol) | (relerr < tol),break;end
end

P=P';