%% Practice problems for exam
%% Problem 1

A = [2 -2 1; 3 2 -3; 1 1 5]; A = A([3 2 1],:)
[L,U,p] = lu(A,'vector')

%
%% Problem 2

L = [1 0 0; -2 1 0; 2 -3 1]; 
U = [2 1 -1; 0 -3 2; 0 0 1];
p = [3; 1; 2];
b = [2; 1; 4];

y = L\b(p)    % solve with forward substitution
x = U\y       % solve with back substitution

%
%% Problem 3

A = [-7 8 -8; -8 9 -9; 6 -7 8];
b = [.97; 1.02; -.01];
xhat = [1; 2; 1];
bhat = A*xhat
epsb = norm(bhat-b,1)/norm(bhat,1)

% Solution 1: We use \frac{\left\|\hat{\;x} -x\right\|}{\left\|\hat{\;x} \right\|}\le \textrm{cond}\left(A\right)\frac{\left\|\hat{\;b} -b\right\|}{\left\|\hat{\;b} \right\|} with \textrm{cond}\left(A\right)=\left\|A\right\|\cdot \left\|A^{-1} \right\|

normA = norm(A,1)
bound = norm(A,1)*25*epsb

% Solution 2: We use A\left(\hat{\;x} -x\right)=\hat{\;b} -b, hence \hat{\;x\;} -x=A^{-1} \left(\hat{\;b} -b\right) and \left\|\hat{\;x\;} -x\right\|\le \left\|A^{-1} \right\|\cdot \left\|\hat{\;b} -b\right\|.
% Therefore \frac{\left\|\hat{\;x\;} -x\right\|}{\left\|\hat{\;x\;} \right\|}\le \frac{\left\|A^{-1} \right\|\cdot \left\|\hat{\;b} -b\right\|}{\left\|\hat{\;x} \right\|}.

bound = 25*norm(bhat-b,1)/norm(xhat,1)

% This is the error bound \frac{\left\|\hat{\;x\;} -x\right\|}{\left\|\hat{\;x\;} \right\|}\le \textrm{cond}\left(A\right)\rho =\frac{\left\|A^{-1} \right\|\cdot \left\|\hat{\;b} -b\right\|}{\left\|\hat{\;x} \right\|}, see "Errors for linear systems".
%
%% Problem 4

t = [0; 1; 3];
y = [2; 4; 3]
A = [t, t.^2]
M = A'*A
v = A'*y
c = M\v    % solve normal equations A*c = v

% Instead of the normal equations we can also use the Matlab \ command:

c = A\y

%
%% Problem 5

f = @(x) x^3+x-11;

% (a) Secant method

x0 = 3; x1 = 2;
x2 = x1 - f(x1)*(x1-x0)/(f(x1)-f(x0))

% (b) Newton method

fp = @(x) 3*x^2 + 1;    % f'(x)
x0 = 2;
x1  = x0 - f(x0)/fp(x0)

% (c) Use fzero: Note that f\left(2\right)=-1 and f\left(3\right)=19 have different signs. Hence there is a solution in the interval \left\lbrack 2,3\right\rbrack. (Since f^{\prime } \left(x\right)=3x^2 +1>0\; the function f is strictly increasing, and the solution is unique.)

xs = fzero(f,[2,3])

%
%% Problem 6
% Define the function f\left(x\right)=\left\lbrack \begin{array}{c}
% f_{1\;} \left(x_1 ,x_2 \right)\\
% f_2 \left(x_{1\;} ,x_2 \right)
% \end{array}\right\rbrack:

f = @(x) [4*x(1)+x(1)^2*x(2)-4; 4*x(2)-x(1)*x(2)^2-3];

% Find the Jacobian f^{\prime } \left(x\right)=\left\lbrack \begin{array}{cc}
% \frac{\partial \;f_1 }{\partial \;x_1 } & \frac{\partial \;f_1 }{\partial \;x_2 }\\
% \frac{\partial \;f_2 }{\partial \;x_1 } & \frac{\partial \;f_2 }{\partial \;x_{2S} }
% \end{array}\right\rbrack

fp = @(x) [4 + 2*x(1)*x(2), x(1)^2; -x(2)^2, 4 - 2*x(1)*x(2)];

% Perform one step of the Newton method starting with x^{\left(0\right)} =\left\lbrack \begin{array}{c}
% 1\\
% 0
% \end{array}\right\rbrack:

x0 = [1;0];
b = f(x0)
A = fp(x0)
d = -A\b        % solve linear system A*d = -b
x1 = x0 + d

% Iterate Newton method until {\left\|f\left(x^{\left(k\right)} \right)\right\|}_1 \le {10}^{-10}:

x = [1;0]
while 1
  b = f(x);
  if norm(b,1)<=1e-10
    break
  end
  A = fp(x);
  d = -A\b;
  x = x + d;
end
x
normf = norm(b,1)

% Use fsolve (with default settings):

xs = fsolve(f,[1;0])
norm(f(xs),1)

%%