Assignment-3-Applied-Math/multi_newton_solver_generalized.m at main · HiGoldAnimationStudios/Assignment-3-Applied-Math · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
%Your Implementation of Newton's method
%INPUTS:
%fun: the mathematical function for which we want to compute the root
% note that the output of fun may include the derivative
% i.e. [fval,dfdx] = fun(x)
% or not, i.e. fval = fun(x)
%x_guess the initial guess for Newton's method
%solver_params: a MATLAB struct containing solver settings
% solver_params.dxtol: terminate early if |x_{i+1}-x_{i}|<dxtol
% solver_params.ftol: terminate early if |f(x_i)|<ftol
% solver_params.dxmax: terminate early if |x_{i+1}-x_{i}|>dxmax
% (meaning that jacobian is probably singular, which is bad)
% solver_params.numerical_diff: boolean (0 or 1)
% 1-> numerically differentiate, 0-> use analytical derivative
%should use numerical differentation to compute the Jacobian or
%use a Jacobian generated by fun. No input assumes numerical by default
% true->fun is assumed to return [fval,J]
% false->fun is assumed to only return fval
%OUTPUTS:
%x: the estimate of the root computed by the function
% exit_flag: an integer indicating whether or not the solver succeeded
function [x, exit_flag, num_evals] = multi_newton_solver_generalized(fun,x_guess,solver_params)
%unpack values from struct (if fields in struct have been set)
    dxtol = 1e-14;
    if isfield(solver_params,'dxtol')
        dxtol = solver_params.dxtol;
    end
    ftol = 1e-14;
    if isfield(solver_params,'ftol')
        ftol = solver_params.ftol;
    end
    dxmax = 1e8;
    if isfield(solver_params,'dxmax')
        dxmax = solver_params.dxmax;
    end
    numerical_diff = 1;
    if isfield(solver_params,'numerical_diff')
        numerical_diff = solver_params.numerical_diff;
    end

    max_iter = 200;
    if isfield(solver_params,'max_iter')
        max_iter = solver_params.max_iter;
    end

    %your code here
    exit_flag=0;

    num_evals=0;

    if numerical_diff==0
        %with J given
        [fval, J] = fun(x_guess);
        num_evals=num_evals+1;
    end

    if numerical_diff==1
        %no J
        fval=fun(x_guess);
        [J,num_evals_1]=approximate_jacobian(fun, x_guess);
        num_evals=num_evals+num_evals_1;
    end

    delta_x = -J\fval;
    count = 0;

    while count<max_iter && norm(delta_x)>dxtol && norm(fval)>ftol && norm(delta_x)<dxmax
        count = count+1;

        if numerical_diff==0
            %with J given
            [fval, J] = fun(x_guess);
            num_evals=num_evals+1;
        end

        if numerical_diff==1
            %no J
            fval=fun(x_guess);
            [J,num_evals_1]=approximate_jacobian(fun, x_guess);
            num_evals=num_evals+num_evals_1;
        end

        delta_x = -J\fval;
        x_guess = x_guess + delta_x;
    end
    % x=newton_Jacobian(fun,J,x_guess);

    if norm(fval)<1e-7
        exit_flag=1;
    end

    x = x_guess;

end

%Implementation of finite difference approximation
%for Jacobian of multidimensional function
%INPUTS:
%fun: the mathetmatical function we want to differentiate
%x: the input value of fun that we want to compute the derivative at
%OUTPUTS:
%J: approximation of Jacobian of fun at x
function [J, num_evals] = approximate_jacobian(fun,x)

    %set jacobian matrix size
    J = zeros(length(fun(x)),length(x));
    num_evals=0;

    for i=1:length(x)
        %calculate column of the jacobian which corresponds to the partial
        %derivatives for given each given set of guesses
        [dfdx,num_evals_i]=approximate_derivative(fun,x,i);
        J(:,i)=dfdx;
        num_evals=num_evals+num_evals_i;
    end
end

%example of how to implement finite difference approximation
%for the first and second derivative of a function
%INPUTS:
%fun: the mathetmatical function we want to differentiate
%x: the input value of fun that we want to compute the derivative at
%OUTPUTS:
%dfdx: approximation of fun'(x)
function [dfdx, num_evals] = approximate_derivative(fun,x,i)

    %set the step size to be tiny
    delta_x = 1e-6;
    ei=zeros(length(x),1);
    ei(i)=1;

    %compute the function at different points near x
    f_left = fun(x-delta_x*ei);
    f_right = fun(x+delta_x*ei);
    num_evals=2;
    %approximate the first derivative
    dfdx = (f_right-f_left)/(2*delta_x);
end