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;