Sorry but cannot fsolve this and don't know why

Good evening MATLAB community,

I am conceding defeat after working at this for the last 3 days and I hope somebody can tell me what I have done incorrectly.

I am trying to extract two values from a set of data, the real and imaginary part of the index of refraction: N=n+ik=n1(1)+i*n1(2). The two equations describing the reflection and transmission coefficients are nonlinear in N so I thought using fsolve would be a good approach. The paper I am working from states that using the Newton method should result in solutions quickly; since fsolve is a variation on the Newton method I thought it was perfectly suited.

There are actually 3 materials in this system where the light goes from air into one medium and then into a second medium. These have indices n0, n1, and n2.

 clear, clc

 n0=[1 0];      %index of refraction for air/first layer
 n10=[1.8 0];%initial guess for index of refraction for second layer
 n2=[3.4 0]; %index of refraction for the third layer

 lambda=10 %wavelength

 d=2.3; %thickness of layer 2

    r=0.297520661157025; %reflection coefficient for Si
    t=0.702479338842975; %transmission coefficient for Si
 next(w,:)=extract(r,t,n0,n2,d, lambda(w),n10); %call for fsolve function

My function extract is

    function y = extract(r, t, n0, n2, d, lambda, n10)

   y = fsolve(@call, n10);

   function y = call(n1) ;

   
  y = [

          -(1+r)...
            +t/(4*n0(1)*n2(1)*(n1(1)^2+n1(2)))*...
            ((n0(1)^2+n1(1)^2+n1(2)^2)*...
            ((n1(1)^2+n2(1)^2+n1(2)^2+n2(2)^2)*cosh(2*2*pi*n1(2)*d/lambda)+...
            2*(n1(1)*n2(1)+n1(2)*n2(2))*sinh(2*2*pi*n1(2)*d/lambda))+...
            (n0(1)^2-n1(1)^2-n1(2)^2)*...
            ((n1(1)^2-n2(1)^2+n1(2)^2-n2(2)^2)*cos(2*2*pi*n1(1)*d/lambda)+...
            2*(n1(1)*n2(2)-n2(1)*n1(2))*sin(2*2*pi*n1(1)*d/lambda))),...

          -(1-r)...
            +t/(2*n2(1)*(n1(1)^2+n1(2)))*...
            (n1(1)*...
            ((n1(1)^2+n2(1)^2+n1(2)^2+n2(2)^2)*sinh(2*2*pi*n1(2)*d/lambda)+...
            2*(n1(1)*n2(1)+n1(2)*n2(2))*cosh(2*2*pi*n1(2)*d/lambda))+...
            n1(2)*...
            ((n1(1)^2-n2(1)^2+n1(2)^2-n2(2)^2)*sin(2*2*pi*n1(1)*d/lambda)-...
            2*(n1(1)*n2(2)-n2(1)*n1(2))*cos(2*2*pi*n1(1)*d/lambda)))

          ];
    end
   end

Now I have verified my equations are correct by solving for r and t for a dimple junction between air and silicon n0=[1 0] and n1=[3.4 0] and n2=[3.4 0]. That is how I obtained r and t in my code above and it follows the simple normal incidence Fresnel coefficients. So my simple test is to insert r and t obtained from the equations already back into the fsolve program to obtain n1=[3.4 0] but it does not work. The warning I get is:

 Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using
 Levenberg-Marquardt algorithm instead. 
 > In fsolve at 324
   In extract at 4
   In myfun_3 at 14
 Optimizer appears to be converging to a minimum that is not a root:
 Sum of squares of the function values exceeds the square root of 
 options.TolFun. Try again with a new starting point.

I have tried to use other peoples forum queries to solve my problem but I have not been able to fix it.

If anybody can take a look at it and give me some insight I would greatly appreciate it and I am sorry for asking something other people have asked about in the past.

Thank you all for your consideration,

Justin

Answer 1

[Sargondjani]

fsolve can only find local solutions. like the exit message says: it finds a local minimum but it is not 0. So you should try a different starting value... Did you try this?

also, it seems that you can calculate the gradient analytically. providing the analytical gradient is always a good idea (it should be the second output argument of the objective function and then you set 'Jacobian' to 'on' in options)

Answer 2

[Teja Muppirala]

Ok, you have a very subtle, but pathological problem with your code.

Compare this:

y = [ (1)...
+2]

To this:

y = [ (1)+...
2]

The first gives you y = [1 2], while the second gives you y = 3. This is exactly the problem in your code. As it is written, y will be a 4-element vector, when you really mean for it to be a 2-element vector. You are essentially trying to solve 4 equations with 2 unknowns, which is why you are getting that message.

This can be fixed by changing this:

-(1+r)...

to this:

-(1+r)+...

in the two places that it occurs. When I do that, FSOLVE converges instantly to

y = 2.173587481230821 -0.000000000000000