\chapter{Soluciones a las actividades}


\section{Actividades del capítulo 2}

\hrulefill

\begin{lstlisting}[language=matlab, title = 2.1.1]
% h^5 debe ser mayor que 10*eps donde el 10 viene del mayor coef.

h = 2*0.0011730;
x = linspace(1-h,1+h,100);
y = x.^5 - 5*x.^4 + 10*x.^3 - 10*x.^2 + 5*x - 1;
plot(x,y)
\end{lstlisting}

\begin{lstlisting}[language=matlab, title = 2.1.2]
% Discretización muy grande. Cuanta mayor es la discretización más son los puntos de impresiónentre los que se interpola y en general dan valores erróneos.

h = 1e-5;
x = linspace(1-h,1+h,1000);
y = x.^5 - 5*x.^4 + 10*x.^3 - 10*x.^2 + 5*x - 1;
plot(x,y)

% Cambio de signo (no hay cancelación)

h = 1e-5;
x = linspace(1-h,1+h,100);
y = x.^5 - 5*x.^4 + 10*x.^3 - 10*x.^2 + 5*x + 1;
plot(x,y)
\end{lstlisting}

\begin{lstlisting}[language=matlab, title = 2.1.3]
% En el primero \sqrt{1+a_{n-1}^2}-1 es aproximadamente
% \sqrt(1+ \pi^2 2^{-2n}) -1 \approx
% (1+\pi^2 2^{-2n-1}) -1
% que debe dar errores relativos grandes cuando 
% \pi^2 2^{-2n-1} \approx eps 
% Despejando con eps = 2^{-52} se tiene n =27.151
\end{lstlisting}

\begin{lstlisting}[language=matlab, title = 2.1.4]
function y = alg1(x)
    % Primera fórmula
    y = (sqrt(1+x^2)-1)/x;
end
\end{lstlisting}


\begin{lstlisting}[language=matlab]
function y = alg2(x)
    % Segunda fórmula
    y = x/(sqrt(1+x^2)+1);
end
\end{lstlisting}

\begin{lstlisting}[language=matlab, title = 2.1.5]
% aproximando tan(x) por x +x^3/3 (Taylor) 
% 2^{n+1}a_n -\pi \approx \pi^3/3/2^{2n+2}

% La predicción funciona bien hasta la iteración 24
anp = 1;
for k = 2:28
    anp = anp/(sqrt(1+anp^2)+1);
    res = 2^(k+1)*anp;
	pred = pi^3/3/2^(2*k+2) ;
    disp(['Iteración ' num2str(k) ' -> ' num2str(res - pi, ' %.10e') ' -> ' num2str(pred, ' %.10e')])
end
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 2.2.1]
% Matriz
A = [1,2,-1; 2, 5, 1; -1,1,11 ];
% Inversa con Gauss-Jordan sobre [A I]
E = rref( [A, eye(3)] );
% Esto muestra las últimas tres columnas
E(:,4:6)

% Inversa con un comando directo (también A^-1)
inv(A)
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 2.3.1]
% Se busca que los pivotes sean grandes
A = [10,1; 1,1];
[L,U] = lu(A);
disp(L)

% Se busca que los pivotes sean grandes
A = [1,10; 3,1];
[L,U] = lu(A);
disp(L)
\end{lstlisting}


\hrulefill


\section{Actividades del capítulo 3}

\hrulefill



\begin{lstlisting}[language=matlab, title = 3.1.1]
% Tomamos por ejemplo
A = [1,2,0;2,4,1;1,0,1]
% En octave el error es 
% warning: division by zero
% En matlab simplifimente aparece  NaN e Inf en las matrices
\end{lstlisting}


\hrulefill

\begin{lstlisting}[language=matlab, title = 3.2.1]
% No funciona con v = k:n
% A(k+1,k) está por debajo de la diagonal y por tanto
% no corresponde a U. Estaríamos sobreescribiendo valores de L.
% En el algortimo original
% L(k+1,k) = U(k+1,k)/U(k,k) (línea 8)
% U(k+1,k) = U(k+1,k) - L(k+1,k)*U(k,k) (línea 10)
% hace que U(k+1,k)=0 y no hay problema. 
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 3.3.1]
function x = ltrs(A, b)
    % Resuelve sistema con matriz triangular inferior
    n = size(A,1);
    x = zeros(n, 1);
    x(1) = b(1)/A(1,1);
    for ii = 2:n
        x(ii) = (b(ii) - A(ii,1:ii-1)*x(1:ii-1) )/A(ii,ii);
    end
end
\end{lstlisting}

\begin{lstlisting}[language=matlab]
function x = utrs(A, b)
    % Resuelve sistema con matriz triangular superior
    n = size(A,1);
    x = zeros(n, 1);
    x(n) = b(n)/A(n,n);
    for ii = n-1:-1:1
        x(ii) = (b(ii) - A(ii,(ii+1):n)*x((ii+1):n) )/A(ii,ii);
    end
end
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 3.3.2]
A = [1,2,1;3,4,5;5,8,1];
b = [4;12;14];
[L,U] = mlu(A);
y = ltrs(L,b);
x = utrs(U,y);
disp(x)
\end{lstlisting}

\begin{lstlisting}[language=matlab]
function [L,U] = mlu(A)
    % Descomposición LU
    n = size(A,1);
    U = A;
    L = eye(n);
    
    for k = 1:n-1
        v = k:n;
        for ii = k+1:n
            L(ii,k) = U(ii,k)/U(k,k);
            U(ii,v) = U(ii,v) - L(ii,k)*U(k,v);
        end
    end
end
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 3.4.1]
N = 1000;
x = linspace(0,4,N);
f = x.^2.*exp(-x);
[m,j] = max(f);
disp(m)
disp(x(j))
\end{lstlisting}


\hrulefill

\begin{lstlisting}[language=matlab, title = 3.4.2]
function [L,U,P] = plu(A)
% Descomposición LU con pivotaje
    n = size(A,1);
    U = A;
    L = eye(n);
    P = eye(n);
    
    for k = 1:n-1
    	[m,j] = max( abs(U(k:n,k)) );
    	
    	if j~=1
    		r = j+k-1;
    		U([k,r], k:n) = U([r,k], k:n);
    		L([k,r], 1:k-1) = L([r,k], 1:k-1);
    		P([k,r], :) = P([r,k], :);
    	end
    	
        v = k:n;
        for ii = k+1:n
            L(ii,k) = U(ii,k)/U(k,k);
            U(ii,v) = U(ii,v) - L(ii,k)*U(k,v);
        end
    end
end
\end{lstlisting}


\hrulefill

\begin{lstlisting}[language=matlab, title = 3.5.1]
% Matlab online
% Elapsed time is 0.099433 seconds.
% Elapsed time is 1.427905 seconds.
%
% Recorrer todas las filas de una columna fijada es más rápido
% => se almacena por columnas, como en Fortran.
\end{lstlisting}

{\tiny\url{https://www.mathworks.com/help/matlab/matlab_external/matlab-data.html}}

\hrulefill

\begin{lstlisting}[language=matlab, title = 3.5.2]
% Matlab online
N = 1000;
A = rand(N);
tic
mlu(A);
toc

tic
mlut(A);
toc
\end{lstlisting}

\begin{lstlisting}[language=matlab]
function [L,U] = mlut(A)
    % Descomposición L^tU^t traspuestas
    n = size(A,1);
    U = transpose(A);
    L = eye(n);
    
    for k = 1:n-1
        v = k:n;
        for ii = k+1:n
            L(k,ii) = U(k,ii)/U(k,k);
            U(v,ii) = U(v,ii) -L(k,ii)*U(v,k);
        end
    end
end
\end{lstlisting}

\begin{lstlisting}[language=matlab]
% Matlab online
% Elapsed time is 7.671028 seconds.
% Elapsed time is 5.264230 seconds.
%
% El algoritmo es más rápido
\end{lstlisting}


\hrulefill



\section{Actividades del capítulo 4}

\hrulefill



\begin{lstlisting}[language=matlab, title = 4.1.1]
N = 10;
A = rand(N) + i*rand(N);
disp( sqrt( max( eig(A'*A) ) ) )
disp( norm(A) )
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 4.1.2]
N = 3;
maxp = 0;
for j= 1:1000
	A = rand(N);
	B = rand(N);
	res = norm(A*B,5)/norm(A,5)/norm(B,5);
	if res>maxp
		maxp = res;
	end
end
disp(['El máximo es ' num2str(maxp)])
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 4.2.1]
format long
p = 1;
q = 1;
for k = 1:6
	ptemp = p^2+2*q^2;
	q = 2*p*q;
	p = ptemp;
	disp(['num/den = ' num2str(p) '/' num2str(q)])
end
% En la sexta iteración ya los números son 
% mayores que 2^{63} y por tanto no fiables.
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 4.2.2]
% Para a<1 converge para todo x_0.
% Para a>1 diverge (para x_0 ~= la solución).
% Para a=1 no hay ecuación.
N = 5;
x = 1;
p1 = [];
p2 = [];
for k = 1:N
	% (x_n,x_{n+1}), (x_{n+1},x_{n+1})
	p1 = [p1,x];
	x = (x+8)/3;
	p1 = [p1,x];
	p2 = [p2,x,x];
end

figure(1)
plot(p1,p2,'--o')
hold on
t = linspace(1,4,10);
plot(t,t, t,(t+8)/3)
hold off
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 4.3.1]
N = 5;
A = [5,2,1; -3,1,1; 2,3,1];
b = [10;6;20];

n = size(A,1);
Ap = A - diag(diag(A));
Dinv = diag(A).^-1;

% Si el radio espectral es menor 
% que 1, no converge en general
r = max(abs(inv(diag(diag(A)))*eig(Ap)))
disp(['Radio espectral = ', num2str(r)])

x = zeros(n,1);
for k = 1:N
	x = Dinv.*(b-Ap*x);
	disp(norm(x))
end
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 4.3.2]
function [x,err] = mjac (A, b, Nmax, tol)
	n = size(A,1);
	Ap = A - diag(diag(A));
	Dinv = diag(A).^-1;

	x = zeros(n,1);
	for k = 1:Nmax
		x_old = x;
		x = Dinv.*(b-Ap*x);
		err = norm(x-x_old)/norm(x_old+eps);
		if err < tol
			return
		end
	end
end
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 4.4.1]
function [x,err] = mgase (A, b, Nmax, tol)
	n = size(A,1);

	x = zeros(n,1);
	for k = 1:Nmax
		x_old = x;
		for ii = 1:n
			x(ii) = (b(ii) - A(ii,1:ii-1)*x(1:ii-1) - A(ii,ii+1:n)*x(ii+1:n) )/A(ii,ii);
		end
		err = norm(x-x_old)/norm(x_old+eps);
		if err < tol
		disp(['Se han usado ',num2str(k),' iteraciones'])
			return
		end
	end
end
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 4.5.1]
N = 10
n = 100
A = rand(n) + (n-1)*eye(n);
b = rand(n,1);
x_exact = A\b;

% xj -> Jacobi
% xgs -> Gauss-Seidel

xgs = zeros(n,1);
xj = zeros(n,1);
Ap = A - diag(diag(A));
Dinv = diag(A).^-1;

for k = 1:N
	xj = Dinv.*(b-Ap*xj);
	for ii = 1:n
		xgs(ii) = (b(ii) - A(ii,1:ii-1)*xgs(1:ii-1) - A(ii,ii+1:n)*xgs(ii+1:n) )/A(ii,ii);
	end
	e1 = norm(xj-x_exact);
	e2 = norm(xgs-x_exact);
	disp( ['k = ' num2str(k) ' -> J: ' num2str(e1,'%.5e') ' -> G-S: ' num2str(e2,'%.5e')]  )
end
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 4.6.1]
% Número iteraciones
N = 100
% Sistema
A = [2,2,1; 8,4,1; 4,7,-4];
b = [6; 14; 3];

% número de filas
m = size(A,1);
% número de comunas
n = size(A,2);

% Vector inicial
x = zeros(n,1);

for k = 1:m
	no = 1/norm(A(k,:));
	A(k,:) = A(k,:)*no;
	b(k) = b(k)*no;
end

for k = 1:N
	ii = mod(k,m) +1;
	fila = A( ii, : );
	x = x + (b(ii) - fila*x)*fila';
end

disp(x)
\end{lstlisting}


\hrulefill


\section{Actividades del capítulo 5}

\hrulefill



\begin{lstlisting}[language=matlab, title = 5.1.1]
N = 300

A = rand(N);

tic
m = size(A,2);
Q = A;
Q(:,1) = A(:,1)/norm(A(:,1));

for k = 2:m
	u_k = A(:,k);
	for j = 1:k-1
		u_k = u_k - Q(:,j)' * A(:,k) * Q(:,j);
	end
	Q(:,k) = u_k/norm( u_k );
end
toc

tic
orth(A);
toc
\end{lstlisting}
Salida (aleatoria) para varios valores de \code{N}:
\begin{lstlisting}[language=matlab, numbers=none]
N =  300
Elapsed time is 0.520701 seconds.
Elapsed time is 0.0822821 seconds.

N =  600
Elapsed time is 2.09508 seconds.
Elapsed time is 0.631554 seconds.

N =  900
Elapsed time is 4.96111 seconds.
Elapsed time is 2.29034 seconds.
\end{lstlisting}


\hrulefill



\begin{lstlisting}[language=matlab, title = 5.1.2]
A = [1,1,0;1,2,0;-1,3,1]
tol = 1e-10;
for j1 = 1:size(A,2)
	for j2 = j1+1:size(A,2)
		if abs(A(:,j1)'*A(:,j2))>tol
			disp([num2str(j1) ' y ' num2str(j2) ' no son ortogonales'])
		else
			disp([num2str(j1) ' y ' num2str(j2) ' son ortogonales'])
		end
	end
end
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 5.2.1]
N = 4
k = 2

v = zeros(N,1);
while v'*v==0
	v = round(k*(2*rand(N,1)-1));
end
d = v'*v;
Q = d*eye(N)-2*v*v';
disp('La matriz es')
disp(Q)
disp(['dividida por ' num2str(d)])
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 5.2.2]
function v = vh(x)
	v = x;
	if v(1)<0
		v(1) = v(1)-norm(x);
	else 
		v(1) = v(1)+norm(x);
	end
	v = v/norm(v);
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 5.3.1]
times1 = [];
times2 = [];
rang = 100:10:500;

for N = rang
	A = rand(N);

	t0 = tic;
	[Q,R] = qr(A);
	times1 = [times1,toc(t0)];
	
	t0 = tic;
	[Q,R] = mQR(A);
	times2 = [times2,toc(t0)];
end

plot( rang, times2./times1)
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 5.3.2]
m = 3
n = 3*m

A = rand(n,m);
[Q,R] = mQR(A);

eig(Q'*Q)
eig(Q*Q')
\end{lstlisting}
Son $m$ unos en el primes caso, completados con $n-m$ ceros en el segundo. 

\hrulefill


\begin{lstlisting}[language=matlab, title = 5.3.3]
times1 = [];
times2 = [];
rang = 100:10:500;

for N = rang
	A = rand(N);

	t0 = tic;
	[Q,R] = qr(A);
	times1 = [times1,toc(t0)];
	
	t0 = tic;
	[Q,R] = mQR2(A);
	times2 = [times2,toc(t0)];
end

plot( rang, times2./times1)
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 5.4.1]
times1 = [];
times2 = [];
rang = 1000:100:3000;

for N = rang
	A = rand(N);
	b = rand(N,1);

	t0 = tic;
	[L,U,P] = lu(A);
	times1 = [times1,toc(t0)];

	t0 = tic;
	[Q,R] = qr(A);
	times2 = [times2,toc(t0)];
end

plot( rang, times2./times1)
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 5.4.2]
N = 10
A = rand(N);
b = rand(N,1);

[L,U,P] = lu(A);
y = ltrs(L,P*b);
x = utrs(U,y);

norm(x-A\b)
\end{lstlisting}

\hrulefill


\section{Actividades del capítulo 6}

\hrulefill


\begin{lstlisting}[language=matlab, title = 6.1.1]
A = [1,2; 1,1; 3,1];
b = [2;8;6];

[Q,R] = qr(A);
m = size(A,2);
T = R(1:m,:);
Qm = Q(:,1:m);

utrs(T, Qm'*b)
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 6.1.2]
function res = rrerr(x,y)
    % Error en la recta de regresión
	n = size(x,1);
	X = [ones(n,1), x];
	c = inv(X'*X)*X'*y;
	res = max( abs(c(1)+c(2)*x-y) );
end
\end{lstlisting}
No se cumple \code{rrerr(x,y)=rrerr(y,x)} en general. 

\begin{lstlisting}[language=matlab]
N = 6
x = (1:N)';
y = 2*x + [1;zeros(N-1,1)];

rrerr(x,y)
rrerr(y,x)
\end{lstlisting}



\hrulefill


\begin{lstlisting}[language=matlab, title = 6.2.1]
N = 10
A = rand(2*N,N)+i*rand(2*N,N);
b = rand(2*N,1)+i*rand(2*N,1);
[U,D,V] = svd(A);
E = diag( diag(D) );
norm( V*inv(E)*U(:,1:N)'*b - A\b )
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 6.2.2]
A = rand(100,200);
res = zeros(99,1);

[U,D,V] = svd(A);

for r = 99:-1:1
	D(r+1,r+1)=0;
	Ar = U*D*V';
	res(r) =  norm(A-Ar,'fro');
end

figure(1)
plot(res);
\end{lstlisting}



\hrulefill



\begin{lstlisting}[language=matlab, title = 6.2.3]
A = [0.546, -0.722; 0.783, 0.427];

x = [2,3,3,2,2; 0,0,1,1,0];
y = A*x;

[U,D,V] = svd(A);

z = U*V'*x;

figure(1)
plot(x(1,:), x(2,:))
hold on
plot(y(1,:), y(2,:))
hold on
plot(z(1,:), z(2,:))
hold off
axis('equal')
\end{lstlisting}


\hrulefill



\begin{lstlisting}[language=matlab, title = 6.2.4]
N = 5;
M = 3;
A = rand(N,M)+i*rand(N,M);

[V,D] = eig(A'*A);
[U,D] = qr(A*V);

norm(A-U*D*V','fro')
\end{lstlisting}



\hrulefill


\section{Actividades del capítulo 7}

\hrulefill

\begin{lstlisting}[language=matlab, title = 7.1.1]
% Crea la matriz
N = 100;
A = diag( ones(N-1,1), 1 );
A = 2*eye(N) + A + A';

% ordena los autovalores
v = sort(eig(A));

% los dibuja
plot(v)
\end{lstlisting}


\hrulefill

\begin{lstlisting}[language=matlab, title = 7.2.1]
A1 = [-7,6,-8,5; -16,10,-13,10; -5,2,-2,3; -8,6,-8,6];
A = -A1;
N = 8;
v = rand(size(A,1),1);
v = v/norm(v);

tol = 0.1;
fl = false;
for k = 1:N
	vn = A*v;
	[m,j] = max(abs(vn));
	v = v/v(j);
	vn = vn/vn(j);
	if norm(v-vn) < tol*norm(vn)
		fl = true;
		v = vn/norm(vn);
		break
	end
	v = vn/norm(vn);
end


if fl
	disp(['Se ha alcanzado la tolerancia exigida con N = ' num2str(k)])
else
	disp('No se ha alcanzado la tolerancia exigida')
end

disp('Aproximación del autovector:')
disp(v)

disp('Aproximación del autovalor:')
disp(v'*A*v)
\end{lstlisting}

\hrulefill



\begin{lstlisting}[language=matlab, title = 7.3.1, numbers=none]
% Basta añadir esto antes de axis('equal')
plot(real(eig(A)),imag(eig(A)),'o');
\end{lstlisting}

\hrulefill



\begin{lstlisting}[language=matlab, title = 7.3.2]
A3 = [5,0,-2,1; 0,5,-1,2; -2,-1,5,0; 1,2,0,5];
R = rand(4);
ep = 0.1;
P = ep*(R+R');
A = A3+P;

[C,D] = eig(A3);
[Cp,D] = eig(A);

disp('Diferencia real de los autovalores')
disp(eig(A)-eig(A3))

disp('Diferencia estimada de los autovalores')
disp( diag(C'*P*C) )
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 7.4.1]
A1 = [-7,6,-8,5; -16,10,-13,10; -5,2,-2,3; -8,6,-8,6];
A = A1;
N = 5;
n = size(A,1);
v = rand(n,1);
v = v/norm(v);
w = v;

la = 0.8+0.8*i;
Ai = inv(A-la*eye(n));

for k = 1:N
	fprintf('Iteración %d\n', k)
	
	% Alg 1
	wn = Ai*w;
	w = wn/norm(wn);
	fprintf('\tError alg1 %.6e\n', abs(w'*A*w-1-i))

	% Alg 2
	vn = (A-la*eye(n))\v;
	v = vn/norm(vn);
	la = v'*A*v;
	fprintf('\tError alg2 %.6e\n', abs(la-1-i))
end
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 7.5.1]
function [C,D] = eigt(A)
    % Autovalores y autovectores de matriz
    % triangular superior.
    % Usa la función utrs
	n = size(A,1);
	C = zeros(n);

	C(1,1) = 1;
	for k=2:n
		C(k,k) = 1;
		T = A(1:k-1, 1:k-1) - A(k,k)*eye(k-1);
		C(1:k-1,k) = utrs( T, -A(1:k-1,k) );
	end

	D = diag(diag(A));
end
\end{lstlisting}

\begin{lstlisting}[language=matlab]
function x = utrs(A, b)
    % Resuelve sistema con matriz triangular superior
    n = size(A,1);
    x = zeros(n, 1);
    x(n) = b(n)/A(n,n);
    for ii = n-1:-1:1
        x(ii) = (b(ii) - A(ii,(ii+1):n)*x((ii+1):n) )/A(ii,ii);
    end
end
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 7.5.2]
% Dimensión
n = 100

% Número de iteraciones
N = 400;

% Matriz simétrica aleatoria
AA = rand(n);
AA = AA+AA';


tic
A = AA;
for k = 1:N
	[Qp,R] = qr(A);
	A =  R*Qp;
end
toc

% error
% disp( norm( sort(diag(A))-sort(eig(AA)) ) )


tic
A = AA;
[P,A] = hess(A);
for k = 1:N
	[Qp,R] = qr(A);
	A =  R*Qp;
end
toc

% error
% disp( norm( sort(diag(A))-sort(eig(AA)) ) )
\end{lstlisting}

\hrulefill

\section{Actividades del capítulo 8}

\hrulefill



\begin{lstlisting}[language=matlab, title = 8.1.1]
N = 5;
c = -6/17
x = 1;
for k = 1:N
	x = c*(x^2 - 2) + x
end
disp(['Error absoluto (con signo) = ' num2str(x-sqrt(2))])
\end{lstlisting}

\hrulefill



\begin{lstlisting}[language=matlab, title = 8.2.1]
function xa = acel (x)
	xa = circshift(x,2) - (circshift(x,1)-circshift(x,2)).^2 ./ (x-2*circshift(x,1)+circshift(x,2));
	xa(1) = x(1);
	xa(2) = x(2);
end
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 8.2.2]
N = 10;

f = @(x) x-(x^2-2)/2;

x = ones(N,1);
x(1) = 1;
a(1) = 1;
for k = 1:N-1
	x(k+1) = f(x(k));
	a(k+1)= a(k)-( f(a(k))-a(k) )^2 / ( f(f(a(k)))-2*f(a(k))+a(k) );
end

for k = 3:N
	disp('----------')
	disp(['Iteración ' num2str(k)])
	xa = x(k-2)- (x(k-1)-x(k-2))^2/(x(k)-2*x(k-1)+x(k-2));
	disp(['Error sin acelerar ' num2str(x(k)-sqrt(2))])
	disp(['Error acelerado ' num2str(xa-sqrt(2))])
	disp(['Error Aitken ' num2str(a(k)-sqrt(2))])
end
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 8.2.3]
N = 4;

f = @(x) (x+2/x)/2;

x = ones(N,1);
x(1) = 1;
a(1) = 1;
for k = 1:N-1
	x(k+1) = f(x(k));
	a(k+1)= a(k)-( f(a(k))-a(k) )^2 / ( f(f(a(k)))-2*f(a(k))+a(k) );
end

for k = 1:N
	disp('----------')
	disp(['Iteración ' num2str(k)])
	disp(['Error algoritmo ' num2str(x(k)-sqrt(2))])
	disp(['Error Aitken ' num2str(a(k)-sqrt(2))])
end
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 8.3.1]
function x = bisec (f,a,b, err)
	while b-a>err
		c = (a+b)/2;
		if f(a)*f(c)<0
			b = c;
		else
			a = c;
		end
	end
	x = (a+b)/2;
end
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 8.4.1]
figure(1)
A = juliak(3, 30, 400, 0.01);
subplot(1,3,1)
imshow(A);

A = juliak(4, 180, 400, 0.1);
subplot(1,3,2)
imshow(A);

subplot(1,3,3)
A = juliak(5, 380, 400, 0.1);
imshow(A);
\end{lstlisting}

\begin{lstlisting}[language=matlab]
function R = juliak (k, Nmax, N, tol)
	% k = degree
	% Nmax = number of steps in Newton-Raphson
	% N = number of pixels
	% tol = tolerance
	
	A = meshgrid(linspace(-3/2,3/2,N),1:N);
	A = A + i*A' + eps;
	z = 2^(1/k)*exp(2*pi*i*[0:k-1]/k);
	
	% Newton-Raphson
	for l = 1:Nmax
		A = (1-1/k)*A+2/k*A.^(1-k);
	end

	R = zeros(N);
	for l = 1:k
		R = R + l*(abs(A-z(l))<tol);
	end
	
	% normalization
	R = im2double(mat2gray(R));
end
\end{lstlisting}

\hrulefill


\begin{lstlisting}[language=matlab, title = 8.5.1]
Nmax = 4;
x = zeros(Nmax,1);
x(1) = 1;
x(2) = 2;
for k = 3:Nmax
	x(k) = x(k-1) - (x(k-1)^2-2)*(x(k-1)-x(k-2)) /(x(k-1)^2-x(k-2)^2);
end

figure(1)
t = linspace(x(1),x(2), 100);
plot( t, t.^2-2, 'g')
hold on 
plot(x, x.^2-2)
hold off
\end{lstlisting}

\hrulefill



\section{Actividades del capítulo 9}

\hrulefill


\begin{lstlisting}[language=matlab, title = 9.1.1]
K = 10
a = -K;
b = K;
x = -K:K;
y = x;
y(K+1) = 0.01;
n = 2*K;
w = zeros(n+1,1);
for j = 0:n
	w(j+1)= (-1)^j*nchoosek(n,j)*2^n/(b-a)^n;
end

figure(1)
plot(x,y,'bo')
hold on
t = linspace(a,b,300);
plot( t,  lagr_bar(x, y, w, t) )
hold off
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 9.1.2]
function w = wj (x)
	n = length(x)-1;
	w = zeros(1,n+1);
	for j = 0:n
		temp = x(j+1)-x;
		temp(j+1) = [];
		w(j+1) = 1/prod(temp);
	end
end
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 9.2.1]
n = 4;
x = linspace(-1,1,n+1);

for j = 0:n
	temp = x;
	temp(j+1) = [];
	Lj = poly(temp)/prod(x(j+1)-temp);

	figure(j+1)
	plot(x,polyval(Lj,x),'bo')
	hold on
	t =  linspace(-1,1,200);
	plot(t, polyval(Lj,t))
	hold off
end
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 9.2.2]
x = -2:2;
y = cos(x)+sin(x);

figure(1)
plot(x,y,'bo')
hold on
t = linspace(-2,2,100);
plot(t,cos(t)+sin(t), 'r--')
hold on
for n = 1:4
	cpi = polyfit(x,y,n);
	plot( t,  polyval(cpi,t) )
	hold on
end
hold off
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 9.3.1]
for n=5:20
	x = linspace(-5,5, n+1);
	y = cos(x);
	cpi = polyfit(x,y,n);

	figure(1)
	t = linspace(-5,5,200);
	plot(t,cos(t), 'r--')
	hold on
	plot( t,  polyval(cpi,t) )
	hold on
end
hold off
\end{lstlisting}


\hrulefill

\begin{lstlisting}[language=matlab, title = 9.3.2]
n = 10;
x = linspace(-5,5, n+1);
g = poly(x);

figure(1)
t = linspace(-5,5,300);
p = abs(polyval(g,t));
plot( t,  p )
\end{lstlisting}


\hrulefill



\begin{lstlisting}[language=matlab, title = 9.4.1]
K = 10;
x = -K:K;
y = x;
y(K+1) = 0.01;

figure(1)
t = linspace(-K,K,200);
sc = spline(x,y);
plot(x,y,'bo', t, ppval(sc,t), '-')
\end{lstlisting}


\hrulefill



\begin{lstlisting}[language=matlab, title = 9.5.1]
J = 4;
P = zeros(6*J+1,2);
P(:,1) = 0:6*J;
P(:,2) = [repmat([0,1,1,0,-1,-1],1,J)'; 0];

figure(1)
n = (size(P,1)-4)/3;
for j = 0:n
	r = 3*j+1:3*j+4;
	plot( P(r,1), P(r,2), 'o', P(r,1), P(r,2), '--')
	hold on
	t = linspace(0,1, 300)';
	R = [(1-t).^3, 3*t.*(1-t).^2, 3*t.^2.*(1-t), t.^3]*P(r,:);
	plot( R(:,1), R(:,2) )
	hold on
end
hold off
\end{lstlisting}


\hrulefill


\begin{lstlisting}[language=matlab, title = 9.5.2]
q1 = @(x) x.^3/6;
q2 = @(x) (x+1).^3/6-2*x.^3/3;
q3 = @(x) q2(1-x);
q4 = @(x) q1(1-x);

t = linspace(0,1, 100);
plot(t-2,q1(t))
hold on
plot(t-1,q2(t))
hold on
plot(t,q3(t))
hold on
plot(t+1,q4(t))

hold off
\end{lstlisting}


\hrulefill