Combining my newton's method function for nonlinear equations and my lu decomposition function.

Hi, I'm trying to combine my newton's method function for nonlinear equations with my LU decomposition function.

This is my function for newton's method:

function x = newtonSysGB(J, fun, x0, xtol, ftol, max_iter)

F_value = fun(x0);

for i = 1:max_iter

delta = J(x0)\-F_value;

Delta_norm = norm(delta);

x0 = x0 + delta;

F_value = fun(x0);

if cond(fun(x0))>10^5

fprintf('Warning!')

end

if cond(fun(x0))>10^15

error('Condition number too high!')

end

F_norm = norm(F_value);

if F_norm < ftol || Delta_norm < xtol

x = x0;

break

end

if F_norm > ftol || Delta_norm > xtol

disp ('Need more iterations!')

end

end

J and fun are two other m-files which are:

function F_vector = fun(x)

F_vector = [x(1)*exp(0.1*x(2)) + 0.1*x(3) - 100; x(1)*exp(0.2*x(2)) + 0.2*x(3) - 120; x(1)*exp(0.3*x(2)) + 0.3*x(3) - 150];

end

function J_matrix = J(x)

J_matrix = zeros(3,3);

J_matrix(1,1) = exp(0.1*x(2)); % df1x1

J_matrix(1,2) = 0.1*exp(0.1*x(2))*x(1); % df1x2

J_matrix(1,3) = 0.1; % df1x3

J_matrix(2,1) = exp(0.2*x(2)); % df2x1

J_matrix(2,2) = 0.2*exp(0.2*x(2))*x(1); % df2x2;

J_matrix(2,3) = 0.2; % df2x3;

J_matrix(3,1) = exp(0.3*x(2)); % df3x1

J_matrix(3,2) = 0.3*exp(0.3*x(2))*x(1); % df3x2

J_matrix(3,3) = 0.3; % df3x3

end

I want to change the bolded part of it, which is just like a system of linear equations to solving it using this LU decomposition function I have.

This is the LU decomposition function that I have:

function X_sol = LUgauss(A,b)

n = length(A);

L = zeros(n);

U = zeros(n);

P = eye(n);

for k=1:n

[~,r] = max(abs(A(k:end,k)));

r = n-(n-k+1)+r;

A([k r],:) = A([r k],:);

P([k r],:) = P([r k],:);

L([k r],:) = L([r k],:);

L(k:n,k) = A(k:n,k) / A(k,k);

U(k,1:n) = A(k,1:n);

A(k+1:n,1:n) = A(k+1:n,1:n) - L(k+1:n,k)*A(k,1:n);

end

U(:,end) = A(:,end);

B = P*b;

n = length(B);

y(1,1) = B(1)/L(1,1);

for i = 2:n

y(i,1) = (B(i)-L(i,1:i-1)*y(1:i-1,1))./L(i,i);

end

n = length(y);

X_sol(n,1) = y(n)/U(n,n);

for i = n-1:-1:1

X_sol(i,1) = (y(i)-U(i,i+1:n)*X_sol(i+1:n,1))./U(i,i);

end

end

THANK YOU SOOO MUCH!

Answer 1

[Fabio Freschi]

Is this what you are asking for?

delta = LUgauss(J(x0),-F_value);

Note that the buil-in lu Matlab function can do the work for you

[Nina]

yes, this is what I was looking for.

I can't use any of the built in functions for this task so I have to make my own lu function.

[Fabio Freschi]

Good! If the answer respont to your original question, please accept it