Problem to solve nonlinear Eq. with Newton's Method

Dear All,

I have problem with convergence using Newtons'method. I am trying to apply Newton's method in Matlab for my problem, and I wrote a code for my problem as following. The goal is to find three variables n, k and d. the expected valus are n=0.288; k=3.536; d=15.83 and the norm of function is 0.04529 and it shows the initial guess is good. But the convergence of calculated data for next itteratiobn is not suitable.

the results presented as below:

Iteration| n | k | d | Error |

0 | 0.20000| 3.50000 | 15.0 | 0.04529 |

Warning: Matrix is close to singular or badly scaled. Results may be inaccurate. RCOND = 8.811893e-18.

> In NEWTON_Nonlinear_nkd2003 (line 89)

1 |-126379826448.59407| -27665191566.24697 | -239395745423.6 | NaN |

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

> In NEWTON_Nonlinear_nkd2003 (line 89)

2 | NaN| NaN | NaN | NaN |

Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND = NaN.

> In NEWTON_Nonlinear_nkd2003 (line 89)

I would appreciate it if anyone could help!

clc; clear; close all;
i=sqrt(-1);
n1=1;

n3=1.515;
lambda=632.8;     %nm
alpha0=-deg2rad([50 55 60]);
teta1=deg2rad(50);

phi0 =deg2rad([75.357 68.897 61.947]);

iteration = 5;

syms n k d;
alpha=alpha0(1);
teta2=asin(n1/(n+i*k)*sin(teta1));
teta3=asin(n1/n3*sin(teta1));
beta=2*pi*d*(n+i*k)*cos(teta2)/lambda;
r12p=(n1*cos(teta2)-(n+i*k)*cos(teta1))/(n1*cos(teta2)+(n+i*k)*cos(teta1));
r23p=((n+i*k)*cos(teta3)-n3*cos(teta2))/((n+i*k)*cos(teta3)+n3*cos(teta2));
r12s=(n1*cos(teta1)-(n+i*k)*cos(teta2))/(n1*cos(teta1)+(n+i*k)*cos(teta2));
r23s=((n+i*k)*cos(teta2)-n3*cos(teta3))/((n+i*k)*cos(teta2)+n3*cos(teta3));
rp=(r12p+r23p*exp(2*i*beta))/(1+r12p*r23p*exp(2*i*beta));
rs=(r12s+r23s*exp(2*i*beta))/(1+r12s*r23s*exp(2*i*beta));
A=0.5*(abs(rp)^2*cos(alpha).^2-abs(rs)^2*sin(alpha).^2);
B=0.5*(rp*conj(rs)+rs*conj(rp))*2*sin(alpha).*cos(alpha);
phi1=atan(B./A)-phi0(1);

alpha=alpha0(2);
teta2=asin(n1/(n+i*k)*sin(teta1));
teta3=asin(n1/n3*sin(teta1));
beta=2*pi*d*(n+i*k)*cos(teta2)/lambda;
r12p=(n1*cos(teta2)-(n+i*k)*cos(teta1))/(n1*cos(teta2)+(n+i*k)*cos(teta1));
r23p=((n+i*k)*cos(teta3)-n3*cos(teta2))/((n+i*k)*cos(teta3)+n3*cos(teta2));
r12s=(n1*cos(teta1)-(n+i*k)*cos(teta2))/(n1*cos(teta1)+(n+i*k)*cos(teta2));
r23s=((n+i*k)*cos(teta2)-n3*cos(teta3))/((n+i*k)*cos(teta2)+n3*cos(teta3));
rp=(r12p+r23p*exp(2*i*beta))/(1+r12p*r23p*exp(2*i*beta));
rs=(r12s+r23s*exp(2*i*beta))/(1+r12s*r23s*exp(2*i*beta));
A=0.5*(abs(rp)^2*cos(alpha).^2-abs(rs)^2*sin(alpha).^2);
B=0.5*(rp*conj(rs)+rs*conj(rp))*2*sin(alpha).*cos(alpha);
phi2=atan(B./A)-phi0(2);

alpha=alpha0(3);
teta2=asin(n1/(n+i*k)*sin(teta1));
teta3=asin(n1/n3*sin(teta1));
beta=2*pi*d*(n+i*k)*cos(teta2)/lambda;
r12p=(n1*cos(teta2)-(n+i*k)*cos(teta1))/(n1*cos(teta2)+(n+i*k)*cos(teta1));
r23p=((n+i*k)*cos(teta3)-n3*cos(teta2))/((n+i*k)*cos(teta3)+n3*cos(teta2));
r12s=(n1*cos(teta1)-(n+i*k)*cos(teta2))/(n1*cos(teta1)+(n+i*k)*cos(teta2));
r23s=((n+i*k)*cos(teta2)-n3*cos(teta3))/((n+i*k)*cos(teta2)+n3*cos(teta3));
rp=(r12p+r23p*exp(2*i*beta))/(1+r12p*r23p*exp(2*i*beta));
rs=(r12s+r23s*exp(2*i*beta))/(1+r12s*r23s*exp(2*i*beta));
A=0.5*(abs(rp)^2*cos(alpha).^2-abs(rs)^2*sin(alpha).^2);
B=0.5*(rp*conj(rs)+rs*conj(rp))*2*sin(alpha).*cos(alpha);
phi3=atan(B./A)-phi0(3);

phi = [phi1 ; phi2 ; phi3];
PHI = matlabFunction(phi,'Vars',{n,k,d});
jacob_phi = [diff(phi1,n) diff(phi1,k) diff(phi1,d); diff(phi2,n) diff(phi2,k) diff(phi2,d); diff(phi3,n) diff(phi3,k) diff(phi3,d) ];
jacob_PHI = matlabFunction(jacob_phi,'Vars',{n,k,d});
error = 10^-5 ;
% Expected result n=0.288;  k=3.536; d=15.83;
n0=0.2; k0=3.5; d0=15;
nkd = [n0 ;k0; d0] ;
PHI1 = feval(PHI, nkd(1), nkd(2),nkd(3));
i = 0;
fprintf('      Iteration|    n     |     k      |     d      |  Error      | \n')
norm1 = norm(PHI1);
fprintf('%10d     |%10.5f| %10.5f | %10.1f | %10.5f |\n',i, n0, k0, d0, norm1)

while true
    
    jacob_PHI1 = feval(jacob_PHI, nkd(1), nkd(2), nkd(3));
    nkd = nkd - jacob_PHI1\PHI1;
    PHI1 = feval(PHI, nkd(1), nkd(2), nkd(3));
    i = i + 1 ;
    norm1 = norm(PHI1);
    fprintf('%10d     |%10.5f| %10.5f | %10.1f | %10.5f |\n',i, nkd(1), nkd(2), nkd(3), norm1)
    
    if norm(PHI1) < error || i == iteration
        break
    end
    
end

[Sam Chak]

Hi @chitua Bella, What are the expected values for ? The expected norm is .

i      = sqrt(-1);
n1     = 1;
n3     = 1.515;
lambda =  632.8;     
alpha0 = -deg2rad([50 55 60]);
teta1  =  deg2rad(50);
phi0   =  deg2rad([75.357 68.897 61.947]);

% iteration = 5;
% syms n k d;

n = 0.288;  
k = 3.536; 
d = 15.83;

alpha = alpha0(1);
teta2 = asin(n1/(n+i*k)*sin(teta1));
teta3 = asin(n1/n3*sin(teta1));
beta  = 2*pi*d*(n+i*k)*cos(teta2)/lambda;
r12p  = (n1*cos(teta2)-(n+i*k)*cos(teta1))/(n1*cos(teta2)+(n+i*k)*cos(teta1));
r23p  = ((n+i*k)*cos(teta3)-n3*cos(teta2))/((n+i*k)*cos(teta3)+n3*cos(teta2));
r12s  = (n1*cos(teta1)-(n+i*k)*cos(teta2))/(n1*cos(teta1)+(n+i*k)*cos(teta2));
r23s  = ((n+i*k)*cos(teta2)-n3*cos(teta3))/((n+i*k)*cos(teta2)+n3*cos(teta3));
rp    = (r12p+r23p*exp(2*i*beta))/(1+r12p*r23p*exp(2*i*beta));
rs    = (r12s+r23s*exp(2*i*beta))/(1+r12s*r23s*exp(2*i*beta));
A     = 0.5*(abs(rp)^2*cos(alpha).^2-abs(rs)^2*sin(alpha).^2);
B     = 0.5*(rp*conj(rs)+rs*conj(rp))*2*sin(alpha).*cos(alpha);
phi1  = atan(B./A)-phi0(1)

phi1 = -0.0174

alpha = alpha0(2);
teta2 = asin(n1/(n+i*k)*sin(teta1));
teta3 = asin(n1/n3*sin(teta1));
beta  = 2*pi*d*(n+i*k)*cos(teta2)/lambda;
r12p  = (n1*cos(teta2)-(n+i*k)*cos(teta1))/(n1*cos(teta2)+(n+i*k)*cos(teta1));
r23p  = ((n+i*k)*cos(teta3)-n3*cos(teta2))/((n+i*k)*cos(teta3)+n3*cos(teta2));
r12s  = (n1*cos(teta1)-(n+i*k)*cos(teta2))/(n1*cos(teta1)+(n+i*k)*cos(teta2));
r23s  = ((n+i*k)*cos(teta2)-n3*cos(teta3))/((n+i*k)*cos(teta2)+n3*cos(teta3));
rp    = (r12p+r23p*exp(2*i*beta))/(1+r12p*r23p*exp(2*i*beta));
rs    = (r12s+r23s*exp(2*i*beta))/(1+r12s*r23s*exp(2*i*beta));
A     = 0.5*(abs(rp)^2*cos(alpha).^2-abs(rs)^2*sin(alpha).^2);
B     = 0.5*(rp*conj(rs)+rs*conj(rp))*2*sin(alpha).*cos(alpha);
phi2  = atan(B./A)-phi0(2)

phi2 = -0.0138

alpha = alpha0(3);
teta2 = asin(n1/(n+i*k)*sin(teta1));
teta3 = asin(n1/n3*sin(teta1));
beta  = 2*pi*d*(n+i*k)*cos(teta2)/lambda;
r12p  = (n1*cos(teta2)-(n+i*k)*cos(teta1))/(n1*cos(teta2)+(n+i*k)*cos(teta1));
r23p  = ((n+i*k)*cos(teta3)-n3*cos(teta2))/((n+i*k)*cos(teta3)+n3*cos(teta2));
r12s  = (n1*cos(teta1)-(n+i*k)*cos(teta2))/(n1*cos(teta1)+(n+i*k)*cos(teta2));
r23s  = ((n+i*k)*cos(teta2)-n3*cos(teta3))/((n+i*k)*cos(teta2)+n3*cos(teta3));
rp    = (r12p+r23p*exp(2*i*beta))/(1+r12p*r23p*exp(2*i*beta));
rs    = (r12s+r23s*exp(2*i*beta))/(1+r12s*r23s*exp(2*i*beta));
A     = 0.5*(abs(rp)^2*cos(alpha).^2-abs(rs)^2*sin(alpha).^2);
B     = 0.5*(rp*conj(rs)+rs*conj(rp))*2*sin(alpha).*cos(alpha);
phi3  = atan(B./A)-phi0(3)

phi3 = -0.0105

phi   = [phi1; phi2; phi3];
norm(phi)   % expected norm is 0.04529

ans = 0.0246

Answer 1

[Torsten]

i=sqrt(-1);
n1=1;
n3=1.515;
lambda=632.8;     %nm
alpha0=-deg2rad([50 55 60]);
teta1=deg2rad(50);
phi0 =deg2rad([75.357 68.897 61.947]);
iteration = 5;
syms n k d;
alpha=alpha0(1);
teta2=asin(n1/(n+i*k)*sin(teta1));
teta3=asin(n1/n3*sin(teta1));
beta=2*pi*d*(n+i*k)*cos(teta2)/lambda;
r12p=(n1*cos(teta2)-(n+i*k)*cos(teta1))/(n1*cos(teta2)+(n+i*k)*cos(teta1));
r23p=((n+i*k)*cos(teta3)-n3*cos(teta2))/((n+i*k)*cos(teta3)+n3*cos(teta2));
r12s=(n1*cos(teta1)-(n+i*k)*cos(teta2))/(n1*cos(teta1)+(n+i*k)*cos(teta2));
r23s=((n+i*k)*cos(teta2)-n3*cos(teta3))/((n+i*k)*cos(teta2)+n3*cos(teta3));
rp=(r12p+r23p*exp(2*i*beta))/(1+r12p*r23p*exp(2*i*beta));
rs=(r12s+r23s*exp(2*i*beta))/(1+r12s*r23s*exp(2*i*beta));
A=0.5*(abs(rp)^2*cos(alpha).^2-abs(rs)^2*sin(alpha).^2);
B=0.5*(rp*conj(rs)+rs*conj(rp))*2*sin(alpha).*cos(alpha);
phi1=atan(B./A)-phi0(1);
alpha=alpha0(2);
teta2=asin(n1/(n+i*k)*sin(teta1));
teta3=asin(n1/n3*sin(teta1));
beta=2*pi*d*(n+i*k)*cos(teta2)/lambda;
r12p=(n1*cos(teta2)-(n+i*k)*cos(teta1))/(n1*cos(teta2)+(n+i*k)*cos(teta1));
r23p=((n+i*k)*cos(teta3)-n3*cos(teta2))/((n+i*k)*cos(teta3)+n3*cos(teta2));
r12s=(n1*cos(teta1)-(n+i*k)*cos(teta2))/(n1*cos(teta1)+(n+i*k)*cos(teta2));
r23s=((n+i*k)*cos(teta2)-n3*cos(teta3))/((n+i*k)*cos(teta2)+n3*cos(teta3));
rp=(r12p+r23p*exp(2*i*beta))/(1+r12p*r23p*exp(2*i*beta));
rs=(r12s+r23s*exp(2*i*beta))/(1+r12s*r23s*exp(2*i*beta));
A=0.5*(abs(rp)^2*cos(alpha).^2-abs(rs)^2*sin(alpha).^2);
B=0.5*(rp*conj(rs)+rs*conj(rp))*2*sin(alpha).*cos(alpha);
phi2=atan(B./A)-phi0(2);
alpha=alpha0(3);
teta2=asin(n1/(n+i*k)*sin(teta1));
teta3=asin(n1/n3*sin(teta1));
beta=2*pi*d*(n+i*k)*cos(teta2)/lambda;
r12p=(n1*cos(teta2)-(n+i*k)*cos(teta1))/(n1*cos(teta2)+(n+i*k)*cos(teta1));
r23p=((n+i*k)*cos(teta3)-n3*cos(teta2))/((n+i*k)*cos(teta3)+n3*cos(teta2));
r12s=(n1*cos(teta1)-(n+i*k)*cos(teta2))/(n1*cos(teta1)+(n+i*k)*cos(teta2));
r23s=((n+i*k)*cos(teta2)-n3*cos(teta3))/((n+i*k)*cos(teta2)+n3*cos(teta3));
rp=(r12p+r23p*exp(2*i*beta))/(1+r12p*r23p*exp(2*i*beta));
rs=(r12s+r23s*exp(2*i*beta))/(1+r12s*r23s*exp(2*i*beta));
A=0.5*(abs(rp)^2*cos(alpha).^2-abs(rs)^2*sin(alpha).^2);
B=0.5*(rp*conj(rs)+rs*conj(rp))*2*sin(alpha).*cos(alpha);
phi3=atan(B./A)-phi0(3);
phi = [phi1 ; phi2 ; phi3];
PHI = matlabFunction(phi,'Vars',{n,k,d});
n0=0.2; k0=3.5; d0=15;
nkd0 = [n0 ;k0; d0] ;
nkd = fsolve(@(x)PHI(x(1),x(2),x(3)),nkd0,optimset('TolFun',1e-8,'TolX',1e-8))

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

nkd = 
  -0.2920 - 0.0471i
   3.5375 + 0.0987i
  15.7937 - 0.7876i

norm(PHI(nkd(1),nkd(2),nkd(3)))

ans = 3.5339e-06