高斯牛顿法：

function [ x_ans ] = GaussNewton( xi, yi, ri)
% input : x = the x vector of 3 points
%         y = the y vector of 3 points
%         r = the radius vector of 3 circles
% output : x_ans = the best answer
    % set up r equations
    r1 = @(x, y) sqrt((x-xi(1))^2+(y-yi(1))^2) - ri(1);
    r2 = @(x, y) sqrt((x-xi(2))^2+(y-yi(2))^2) - ri(2);
    r3 = @(x, y) sqrt((x-xi(3))^2+(y-yi(3))^2) - ri(3);
    % set up Dr matrix
    Dr = @(x) [(x(1) - xi(1))/(sqrt((x(1) - xi(1))^2+(x(2)-yi(1))^2)), (x(2) - yi(1))/(sqrt((x(1) - xi(1))^2+(x(2)-yi(1))^2));
               (x(1) - xi(2))/(sqrt((x(1) - xi(2))^2+(x(2)-yi(2))^2)), (x(2) - yi(2))/(sqrt((x(1) - xi(2))^2+(x(2)-yi(2))^2));
               (x(1) - xi(3))/(sqrt((x(1) - xi(3))^2+(x(2)-yi(3))^2)), (x(2) - yi(3))/(sqrt((x(1) - xi(3))^2+(x(2)-yi(3))^2))];
    % set up r matrix
    r = @(x) [r1(x(1), x(2)); r2(x(1), x(2)); r3(x(1), x(2))];
    x0 = [0, 0]; % initial guess
    while 1
       A = Dr(x0);
       v0 = (A' * A) \ (- A' * r(x0));
       x1 = x0 + v0';
       if norm(x1-x0)<1e-6 % break squest
           break;
       end
       x0 = x1;
    end
    x_ans = x1;
end

Levenberg–Marquardt方法：

function [ x_ans ] = LeveMarq( ti, yi, x_guess, lmd)
% input : t = the x vector of 5 points
%         y = the y vector of 5 points
%         x_guess = the guess vector of x_ans
% output : x_ans = the best answer
    % set up r matrix
    r = @(x) [x(1) * exp(-x(2)*(ti(1) - x(3))^2) - yi(1);
              x(1) * exp(-x(2)*(ti(2) - x(3))^2) - yi(2);
              x(1) * exp(-x(2)*(ti(3) - x(3))^2) - yi(3);
              x(1) * exp(-x(2)*(ti(4) - x(3))^2) - yi(4);
              x(1) * exp(-x(2)*(ti(5) - x(3))^2) - yi(5)];
    % set up Dr matrix
    Dr = @(x) [exp(-x(2)*(ti(1)-x(3))^2), -x(1)*(ti(1)-x(3))^2*exp(-x(2)*(ti(1)-x(3))^2), 2*x(1)*x(2)*(ti(1)-x(3))*exp(-x(2)*(ti(1)-x(3))^2);
               exp(-x(2)*(ti(2)-x(3))^2), -x(1)*(ti(2)-x(3))^2*exp(-x(2)*(ti(2)-x(3))^2), 2*x(1)*x(2)*(ti(2)-x(3))*exp(-x(2)*(ti(2)-x(3))^2);
               exp(-x(2)*(ti(3)-x(3))^2), -x(1)*(ti(3)-x(3))^2*exp(-x(2)*(ti(3)-x(3))^2), 2*x(1)*x(2)*(ti(3)-x(3))*exp(-x(2)*(ti(3)-x(3))^2);
               exp(-x(2)*(ti(4)-x(3))^2), -x(1)*(ti(4)-x(3))^2*exp(-x(2)*(ti(4)-x(3))^2), 2*x(1)*x(2)*(ti(4)-x(3))*exp(-x(2)*(ti(4)-x(3))^2);
               exp(-x(2)*(ti(5)-x(3))^2), -x(1)*(ti(5)-x(3))^2*exp(-x(2)*(ti(5)-x(3))^2), 2*x(1)*x(2)*(ti(5)-x(3))*exp(-x(2)*(ti(5)-x(3))^2)];

    x0 = x_guess; % initial guess
    while 1
       A = Dr(x0);
       M_A = A' * A + lmd .* diag(diag(A' * A));
       M_b = - A' * r(x0);
       v0 = M_A \ M_b;
       x1 = x0 + v0';
       if norm(x1-x0)<1e-6 % break squest
           break;
       end
       x0 = x1;
    end
    x_ans = x1;
end