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

5 views (last 30 days)
Gigi on 8 Nov 2019
Commented: Fabio Freschi on 8 Nov 2019
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!

Fabio Freschi on 8 Nov 2019
Edited: Fabio Freschi on 8 Nov 2019
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