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

Question

Alwaleed Alhazmi 2021-6-28

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/866975-what-is-the-problem-in-my-code-lambert-w-function

回答： Himanshu 2024-6-4

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:

https://www.sciencedirect.com/science/article/pii/S0038092X14006094#t0010

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Himanshu 2024-6-4

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/866975-what-is-the-problem-in-my-code-lambert-w-function#answer_1467136

在 MATLAB Online 中打开

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);