# 數值分析 期中
## Solving nonlinear eqn
### Bisection

%% bisection method
% [root,step] = SolveRootBisection(f,X,type,var)
% f is the interested function.
% X is the initial condition, an array consists of 2 values that bound the 
%     lower and upper limits of the interested root.
% type is the string that determines how the iteration works. There are 2
%     valid inputs, 'TOL' or 'step'.
% var is also determines how the iteration works. When type = 'TOL', 
%     iteration goes on until tolerence < var; if type = 'step', iteration 
%     repeats exactly var times.
% root is the root found.
% step is the step used to reach the tolerence. When type = 'step', step =
%     var.
%%
function [root,step] = SolveRootBisection(f,X,type,var)
    x1 = min(X); x2 = max(X);
    if double(f(x1))*double(f(x2)) < 0 % determine whether the input is valid
        if strcmp(type,'TOL')
            TOL  = var;
            step = 0;
            while abs(x2-x1) > TOL
                step = step+1;
                x3 = (x1+x2)/2;
                if double(f(x1))*double(f(x3)) < 0
                    x2 = x3;
                else
                    x1 = x3;
                end
            end
            root = (x1+x2)/2;
        elseif strcmp(type,'step')
            step = floor(var);
            for r = 1:step
                x3 = (x1+x2)/2;
                if double(f(x1))*double(f(x3)) < 0
                    x2 = x3;
                else
                    x1 = x3;
                end
            end
            root = (x1+x2)/2;
        else
            error('Invalid input');
        end
    else
        error('Invalid input');
    end
end

### Newton method

%% newton method
% [root,step] = SolveRootNewton(f,x1,type,var)
% f is the interested function.
% x is the initial condition.
% type is the string that determines how the iteration works. There are 2
%     valid inputs, 'TOL' or 'step'.
% var is also determines how the iteration works. When type = 'TOL', 
%     iteration goes on until tolerence < var; if type = 'step', iteration 
%     repeats exactly var times.
% root is the root found.
% step is the step used to reach the tolerence. When type = 'step', step =
%     var.
%%
function [root,step] = SolveRootNewton(f,x1,type,var)
    syms f1(x)
        f1(x) = diff(f(x));
    x2 = x1-double(f(x1))/double(f1(x1));
    if strcmp(type,'TOL')
        TOL  = var;
        step = 1;
        while abs(x2-x1) > TOL
            step = step+1;
            x1 = x2;
            x2 = x1-double(f(x1))/double(f1(x1));
        end
    elseif strcmp(type,'step')
        step = var;
        for r = 2:step
            x1 = x2;
            x2 = x1-double(f(x1))/double(f1(x1));
        end
    else
        error('Invalid input');
    end
    root = x2;
end

### Secant method

%% secant method
% [bisection_root,step] = SolveRootSecant(f,X,type,var)
% f is the interested function.
% X is the initial condition, an array consists of 2 values.
% type is the string that determines how the iteration works. There are 2
%     valid inputs, 'TOL' or 'step'.
% var is also determines how the iteration works. When type = 'TOL', 
%     iteration goes on until tolerence < var; if type = 'step', iteration 
%     repeats exactly var times.
% root is the root found.
% step is the step used to reach the tolerence. When type = 'step', step =
%     var.
%%
function [root,step] = SolveRootSecant(f,X,type,var)
    x0 = X(1); x1 = X(2);
    x2 = 0;
    if strcmp(type,'TOL')
        TOL  = var;
        step = 0;
        while abs(x1-x0) > TOL || step == 0
            step = step+1;
            x2 = x1-double(f(x1))*(x0-x1)/(double(f(x0))-double(f(x1)));
            x0 = x1;
            x1 = x2;
        end
    elseif strcmp(type,'step')
        step = var;
        for r = 1:step
            x2 = x1-double(f(x1))*(x0-x1)/(double(f(x0))-double(f(x1)));
            x0 = x1;
            x1 = x2;
        end
    else
        error('Invalid input');
    end
    root = x2;
end

### All in one

%% Root finding with different methods
% [root,step] = [root,step] = SolveRoot(f,method,X,type,var)
% f is the interested function.
% method is the method the user wants to use to find the root
% X is the initial condition, an array consists of 2 values that bound the 
%     lower and upper limits of the interested root.
% type is the string that determines how the iteration works. There are 2
%     valid inputs, 'TOL' or 'step'.
% var is also determines how the iteration works. When type = 'TOL', 
%     iteration goes on until tolerence < var; if type = 'step', iteration 
%     repeats exactly var times.
% root is the root found.
% step is the step used to reach the tolerence. When type = 'step', step =
%     var.
%%
function [root,step] = SolveRoot(f,method,X,type,var)
    if strcmp(method,'bisection')
        x1 = min(X); x2 = max(X);
        if double(f(x1))*double(f(x2)) < 0 % determine whether the input is valid
            if strcmp(type,'TOL')
                TOL  = var;
                step = 0;
                dig  = -floor(log10(TOL));
                while abs(x2-x1) > TOL
                    step = step+1;
                    x3 = (x1+x2)/2;
                    if double(f(x1))*double(f(x3)) < 0
                        x2 = x3;
                    else
                        x1 = x3;
                    end
                end
                root = round((x1+x2)/2,dig);
            elseif strcmp(type,'step')
                step = floor(var);
                for r = 1:step
                    x3 = (x1+x2)/2;
                    if double(f(x1))*double(f(x3)) < 0
                        x2 = x3;
                    else
                        x1 = x3;
                    end
                end
                root = (x1+x2)/2;
            else
                error('Invalid input');
            end
        else
            error('Invalid input');
        end
    elseif strcmp(method,'newton')
        syms f1(x)
        f1(x) = diff(f(x));
        x2 = X-double(f(X))/double(f1(X));
        if strcmp(type,'TOL')
            TOL  = var;
            step = 1;
            dig  = -floor(log10(TOL));
            while abs(x2-X) > TOL
                step = step+1;
                X = x2;
                x2 = X-double(f(X))/double(f1(X));
            end
            root = round(x2,dig);
        elseif strcmp(type,'step')
            step = var;
            for r = 2:step
                X = x2;
                x2 = X-double(f(X))/double(f1(X));
            end
            root = x2;
        else
            error('Invalid input');
        end
    elseif strcmp(method,'secant')
        x0 = X(1); x1 = X(2);
        x2 = 0;
        if strcmp(type,'TOL')
            TOL  = var;
            step = 0;
            dig  = -floor(log10(TOL));
            while abs(x1-x0) > TOL || step == 0
                step = step+1;
                x2 = x1-double(f(x1))*(x0-x1)/(double(f(x0))-double(f(x1)));
                x0 = x1;
                x1 = x2;
            end
            root = round(x2,dig);
        elseif strcmp(type,'step')
            step = var;
            for r = 1:step
                x2 = x1-double(f(x1))*(x0-x1)/(double(f(x0))-double(f(x1)));
                x0 = x1;
                x1 = x2;
            end
            root = x2;
        else
            error('Invalid input');
        end
    end
end

### Nonlinear eqn system (Newton method)

%% Solution of a system of nonlinear equations using Newton's method
% [root,step] = SolveEqnSet(f1,f2,x0,y0,type,var)
% Symbolic Math Toolbox™ is required
% f1 & f2 are the interested functions.
% X is the initial condition, an array consists of 2 values.
% type is the string that determines how the iteration works. There are 2
%     valid inputs, 'TOL' or 'step'.
% var is also determines how the iteration works. When type = 'TOL', 
%     iteration goes on until tolerence < var; if type = 'step', iteration 
%     repeats exactly var times.
% sol is the solution found.
% step is the step used to reach the tolerence. When type = 'step', step =
%     var.
%%
function [sol,step] = SolveEqnSys(f1,f2,X,type,var)
    syms x y
    x0 = X(1); y0 = X(2);
    x2 = x0; y2 = y0; x1 = 0; y1 = 0;
    if strcmp(type,'TOL')
        TOL  = var;
        step = 0;
        while sqrt((x1-x2)^2+(y1-y2)^2) > TOL || step == 0
            step = step+1;
            x1 = x2; y1 = y2;
            JF = det(jacobian([f1,f2],[x,y]));
            JF = double(subs(JF,{x,y},{x1,y1}));
            deltax = (-f1(x1,y1)*diff(f2,y)+f2(x1,y1)*diff(f1,y))/JF;
            deltax = double(subs(deltax,{x,y},{x1,y1}));
            deltay = (-f2(x1,y1)*diff(f1,x)+f1(x1,y1)*diff(f1,y))/JF;
            deltay = double(subs(deltay,{x,y},{x1,y1}));
            x2 = x1+deltax; y2 = y1+deltay; 
        end
    elseif strcmp(type,'step')
        step = var;
        for r = 1:step
            x1 = x2; y1 = y2;
            JF = det(jacobian([f1,f2],[x,y]));
            JF = double(subs(JF,{x,y},{x1,y1}));
            deltax = (-f1(x1,y1)*diff(f2,y)+f2(x1,y1)*diff(f1,y))/JF;
            deltax = double(subs(deltax,{x,y},{x1,y1}));
            deltay = (-f2(x1,y1)*diff(f1,x)+f1(x1,y1)*diff(f1,y))/JF;
            deltay = double(subs(deltay,{x,y},{x1,y1}));
            x2 = x1+deltax; y2 = y1+deltay; 
        end
    end
    sol = [x2,y2];
end

## Solving linear eqn set
### Guass elimination

%% Solving linear equation set with Gauss elimination
function X = SolveLEqnSet_GE(A,B)
    n = size(A,1);
    X = zeros(n,1);
    A = [A B]; % form the augmented matrix
    for i = 1:n-1
        for k = i:n
            if A(k,i) ~= 0
                break
            end
        end
        S = A(k,:);
        A(k,:) = A(i,:);
        A(i,:) = S;
        for j = i+1:n
            A(j,:) = A(j,:)-A(i,:)*(A(j,i)/A(i,i));
        end
    end
    for i = n:-1:1
        M = A(i,n+1);
        for j = i+1:n
            M = M-X(j)*A(i,j);
        end
        X(i) = M/A(i,i); 
    end
end

### Gauss Jordan

%% Solving linear equation set with Gauss-Jordan elimination
function X = SolveLEqnSet_GJ(A,B)
    n = size(A,1);
    A = [A B];
    for i = 1:n
        for k = i:n
            if A(k,i) ~= 0
                break
            end
        end
        S = A(k,:);
        A(k,:) = A(i,:);
        A(i,:) = S;
        A(i,:) = A(i,:)/A(i,i);
        for j = 1:n
            if j ~= i
                A(j,:) = A(j,:)-A(i,:)*(A(j,i)/A(i,i));
            end
        end
    end
    X = A(:,end);
end

### Jacobi method

%% Solving linear equation set with Jacobi iterative method
function [X,step,E] = SolveLEqnSet_J(A,B,type,var)
    X    = zeros(length(B),1); % the solution array
    E    = ones(length(B),1); % the error array
    if strcmp(type,'TOL')
        TOL  = var;
        step = 0;
        while (max(E)>TOL)
            step = step+1;
            Xm   = X;
            for i = 1:length(B)
                M = 0;
                for j = 1:length(B)
                    if j ~= i
                        M = M+A(i,j)*Xm(j);
                    end
                end
                X(i) = (B(i)-M)/A(i,i);
            end
            E = abs((X-Xm)./X);
        end
    elseif strcmp(type,'step')
        step = var;
        for r = 1:step
            Xm = X;
            for i = 1:length(B)
                M = 0;
                for j = 1:length(B)
                    if j ~= i
                        M = M+A(i,j)*Xm(j);
                    end
                end
                X(i) = (B(i)-M)/A(i,i);
            end
            E = abs((X-Xm)./X);
        end
    else
        error('Invalid input');
    end
end

### Guass-Seidel method

%% Solving linear equation set with Gauss-Seidel iterative method
function [X,step,E] = SolveLEqnSet_GS(A,B,type,var)
    X    = zeros(length(B),1); % the solution array
    E    = ones(length(B),1); % the error array
    if strcmp(type,'TOL')
        TOL  = var;
        step = 0;
        while (max(E)>TOL)
            step = step+1;
            for i = 1:length(B)
                Xm = X(i);
                M  = 0;
                for j = 1:length(B)
                    if j ~= i
                        M = M+A(i,j)*X(j);
                    end
                end
                X(i) = (B(i)-M)/A(i,i);
                E(i) = abs((X(i)-Xm)/X(i));
            end
        end
    elseif strcmp(type,'step')
        step = var;
        for r = 1:step
            for i = 1:length(B)
                Xm = X(i);
                M  = 0;
                for j = 1:length(B)
                    if j ~= i
                        M = M+A(i,j)*X(j);
                    end
                end
                X(i) = (B(i)-M)/A(i,i);
                E(i) = abs((X(i)-Xm)/X(i));
            end
        end
    else
        error('Invalid input');
    end
end

## Interpolation
### Linear regression

%% Linear regression
function [C,f] = PolyFit1(x,y)
    n   = length(x);
    x   = reshape(x,1,n); % ensure products function as expected
    y   = reshape(y,1,n);
    Sx  = sum(x);
    Sy  = sum(y);
    Sxy = sum(x.*y);
    Sxx = sum(x.*x);
    C0  = (Sxx*Sy-Sxy*Sx)/(n*Sxx-(Sx)^2);
    C1  = (n*Sxy-Sx*Sy)/(n*Sxx-(Sx)^2);
    C   = [C0,C1];
    syms f(x)
    f(x) = C(1)+x*C(2);
end

### Polynomial regression

%% mth-order polynomial regression
function [C,f] = PolyFitm(x,y,m)
    n = length(x);
    x = reshape(x,1,n); % ensure products function as expected
    y = reshape(y,1,n);
    A = zeros(m+1,m+1);
    B = zeros(m+1,1);
    for i = 1:m+1
        for j = 1:m+1
            A(i,j) = sum(x.^(i+j-2));
        end
        B(i,1) = sum(x.^(i-1).*y);
    end
    C = SolveLEqnSet_GJ(A,B);
    syms f(x)
    f(x) = 0;
    for n = 1:m+1
        f(x) = f(x)+C(n)*x.^(n-1);
    end
end