# What is the problem in my code? Lambert-W function

4 views (last 30 days)
Alwaleed Alhazmi on 28 Jun 2021
Answered: Himanshu on 4 Jun 2024
I'm trying to solve a non-linear system of equations using Newton-Raphson method. The problem is in section 2.1 from this paper:
It goes only up to the second iteration, and I got incorrect answer.
I used these values for iN and VN (N=1,2,3,4,5). V1 = 1.00400400e+00; V2 = 1.02598500e+00; V3 = 1.05325500e+00; V4 = 1.08812300e+00; V5 = 1.13388700e+00; i1 = -2.40036700e-02; i2 = -3.59849700e-02; i3 = -5.32552100e-02; i4 = -7.81225600e-02 ; i5 = -1.13887200e-01;
I used a = Iph, b = I0, f = n, c = Rs, d = Rsh, Vth = 0.0258
Here's my code:
clear all; clc;
N = 20;
tErr = 1e-6;
syms a b c d f
V1 = 1.00400400e+00; V2 = 1.02598500e+00; V3 = 1.05325500e+00; V4 = 1.08812300e+00; V5 = 1.13388700e+00; i1 = -2.40036700e-02; i2 = -3.59849700e-02; i3 = -5.32552100e-02; i4 = -7.81225600e-02 ; i5 = -1.13887200e-01;
m1 = d*(-i1 + a + b)/(0.0258*f); m2 = d*(-i2 + a + b)/(0.0258*f); m3 = d*(-i3 + a + b)/(0.0258*f); m4 = d*(-i4 + a + b)/(0.0258*f);
m5 = d*(-i5 + a + b)/(0.0258*f); n = .0258*f; k = b*d/n;
f1 = -i1*c-i1*d+a*d-n*lambertw(b*d*exp(m1)/n)-V1;
f2 = -i2*c-i2*d+a*d-n*lambertw(b*d*exp(m2)/n)-V2;
f3 = -i3*c-i3*d+a*d-n*lambertw(b*d*exp(m3)/n)-V3;
f4 = -i4*c-i4*d+a*d-n*lambertw(b*d*exp(m4)/n)-V4;
f5 = -i5*c-i5*d+a*d-n*lambertw(b*d*exp(m5)/n)-V5;
fs = [f1; f2; f3; f4; f5];
xs = [a; b; f; c; d];
jac = [d/lambertw(b*d*exp(m1)/n), -(n*lambertw(b*d*exp(m1)/n)-b*d)/(b*(1+lambertw(b*d*exp(m1)/n))), lambertw(b*d*exp(m1)/n)*(-n*lambertw(b*d*exp(m1)/n)-i1*d+a*d+b*d)/(f*(1+lambertw(b*d*exp(m1)/n))), -i1, -(n*lambertw(b*d*exp(m1)/n)+i1*d-a*d-b*d)/((1+lambertw(b*d*exp(m1)/n))*d);
d/lambertw(b*d*exp(m2)/n), -(n*lambertw(b*d*exp(m2)/n)-b*d)/(b*(1+lambertw(b*d*exp(m2)/n))), lambertw(b*d*exp(m2)/n)*(-n*lambertw(b*d*exp(m2)/n)-i2*d+a*d+b*d)/(f*(1+lambertw(b*d*exp(m2)/n))), -i2, -(n*lambertw(b*d*exp(m2)/n)+i2*d-a*d-b*d)/((1+lambertw(b*d*exp(m2)/n))*d);
d/lambertw(b*d*exp(m3)/n), -(n*lambertw(b*d*exp(m3)/n)-b*d)/(b*(1+lambertw(b*d*exp(m3)/n))), lambertw(b*d*exp(m3)/n)*(-n*lambertw(b*d*exp(m3)/n)-i3*d+a*d+b*d)/(f*(1+lambertw(b*d*exp(m3)/n))), -i3, -(n*lambertw(b*d*exp(m3)/n)+i3*d-a*d-b*d)/((1+lambertw(b*d*exp(m3)/n))*d);
d/lambertw(b*d*exp(m4)/n), -(n*lambertw(b*d*exp(m4)/n)-b*d)/(b*(1+lambertw(b*d*exp(m4)/n))), lambertw(b*d*exp(m4)/n)*(-n*lambertw(b*d*exp(m4)/n)-i4*d+a*d+b*d)/(f*(1+lambertw(b*d*exp(m4)/n))), -i4, -(n*lambertw(b*d*exp(m4)/n)+i4*d-a*d-b*d)/((1+lambertw(b*d*exp(m4)/n))*d);
d/lambertw(b*d*exp(m5)/n), -(n*lambertw(b*d*exp(m5)/n)-b*d)/(b*(1+lambertw(b*d*exp(m5)/n))), lambertw(b*d*exp(m5)/n)*(-n*lambertw(b*d*exp(m5)/n)-i5*d+a*d+b*d)/(f*(1+lambertw(b*d*exp(m5)/n))), -i5, -(n*lambertw(b*d*exp(m5)/n)+i5*d-a*d-b*d)/((1+lambertw(b*d*exp(m5)/n))*d)];
x = [1.2; 0.001; 1.5; 1; 15];
for i=1:N
f = double(subs(fs, xs, x));
j = double(subs(jac, xs, x));
x = x - inv(j)*f;
err = abs(inv(j)*f);
if(err<tErr)
break;
end
end
display(x);
##### 0 CommentsShow -2 older commentsHide -2 older comments

Sign in to comment.

### Answers (1)

Himanshu on 4 Jun 2024
Hey,
Your implementation of the Newton-Raphson method seems correct in its approach, but there are a few potential changes that could lead to the correct answer:
• Inversion of Jacobian: Directly using inv(j) to compute the inverse of the Jacobian matrix can lead to numerical instability or inaccuracies, especially if the Jacobian is close to singular. It's generally better to solve the linear system (J \Delta x = -f) using more stable methods like x = x - j\f; which internally uses matrix decomposition methods for solving linear equations.
• Error Calculation: The way you calculate the error err = abs(inv(j)*f); is not typical for the Newton-Raphson method. Normally, the error is calculated as the norm of the difference between successive iterations or the norm of the function values. For example, using norm(f) or norm(x - x_prev) where x_prev is the value of x from the previous iteration.
Please refer to the followinf code, which is a slightly modified version of you loop keeing in mind the above mentioned points.
for i=1:N
f_val = double(subs(fs, xs, x));
j_val = double(subs(jac, xs, x));
dx = -j_val\f_val; % More stable way to solve J dx = -f
x = x + dx;
err = norm(dx); % Using norm of dx as the convergence criterion
if(err < tErr)
break;
end
end
display(x);
##### 0 CommentsShow -2 older commentsHide -2 older comments

Sign in to comment.

R2021a

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!