Solving system of nonlinear equations
proland
Mon, 04/04/2011 - 08:51 am
I am currently trying to solve a system of nonlinear equations (2 equations, 2 unknowns). I am trying to use Newton's method, but if fails when the slope is small (an inherent flaw in many iterative methods for solving systems of nonlinear equations). The slope is currently a linear approximation due to the complexity of the equations. Any advice on other methods i could attempt to solve a system of nonlinear equations?
Additional info:
The equations I'm trying to solve, see attached image, are for transmission and reflection of light from a thin film on glass. The variables to solve for, nf and kf, are terms of the complex index of refraction for the particular wavelength being evaluated.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
April 4, 2011 at 09:57 am - Permalink
The function root() uses the findroots function for 2 nonlinear equations, the macro NewtonRaphson calculates the values based on an iterative process. The saved file compares these results in graph7.
While the findroots function is better, it still has a point where it fails. The v_flag value is 1, "user abort", even though i was given the error "rootfinder is not making progress".
I've tried changed the value of /f (trust region factor), presenting the best results in the attached file. Sory for the unorganized procedure window, I've been trying various solutions.
Any thoughts on how to proceed?
Thank you.
April 4, 2011 at 12:34 pm - Permalink
Hm. That sounds like a bug. I'd like to reproduce that.
In order to try it out, I need to be able to run some stuff. Can you give me a tour of the code?
One thing I'd like to do is to make a contour plot of function_f and function_g in the area around a presumed root. Can you tell me how to do that?
I see that there is a loop that runs over 418 iterations, and one iteration uses the previous solution as the starting guess. I also see a couple graphs with a sharp discontinuities. Is it possible that you're trying to use a solution from one side as the starting guess for the other side?
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
April 5, 2011 at 09:45 am - Permalink
The 418 iterations are to calculate n, k at the 418 wavelengths for which I measured reflection and transmission. I am using the previous solution as a first guess for the current iteration becuase the equations for n, k should be continuous.
I've written a short function do make the contour plots you've asked for, function_f and function_g around the area of the supposed root (400 possible values for n, k centered around the guessed root are used to construct a 400 x 400 matrix). These contour plots are the calculated values of the equations used in the two functions.
Just use the function contour(z) where z is the data point to look at. While doing this, I noticed that the contour plots become quite chaotic when looking at data points above 355. Probably a good indication of why my method encountered a discontinuity here.
April 5, 2011 at 01:20 pm - Permalink
It not only gets chaotic, it goes through a region where there isn't any solution, at least within the range of parameters that you're looking at. I made a contour plot using your two matrices fcontour and gcontour, but put them on the same graph. I colored the zero contour black so that it stands out. The plot made by contour(356) is attached. You can see there is no intersection between the two zero contours.
The situation around contour(400) is OK- there is a discontinuity in the surfaces that is probably far enough away from the solution that with the right starting guess it would probably succeed. But the starting guess would pretty much have to be on the correct side of the discontinuity.
And the picture at contour(358) is just total chaos. You could probably get your solutions excluding that chaotic range by working from either end of the range toward the chaotic middle.
I added a function to your experiment file, DoManyContours(start, end), that steps through a range of values so you can watch a little movie of the changes in the contours. It is instructive to watch it pass through chaos and then recover :)
I have attached my experiment file as well as the PNG image. I had to zip it because Igor Exchange thought the unzipped file was a bit too big.
John Weeks
WaveMetrics, Inc.
support@wavemetrics.com
April 5, 2011 at 03:26 pm - Permalink
Seriously though, I do exactly the same kind of work, but with X-rays and neutrons. We measure reflectivity (Amplitude squared) as a function of momentum transfer (depends on wavelength and angle of incidence). We analyse the data with the Abeles transfer matrix method (see the Abeles XOP on IGORexchange). That XOP won't be entirely applicable to you. However, you should be able to work out the reflectivity and transmissivity of the system by simply using the fresnel equations - they are shown in the thumbnail image you provided. Multilayers add a short amount of complexity.
this means you have to provide inputs for each data point (assuming your system has no roughness): wavelength of light, angle of incidence, and refractive indices and thicknesses of the fronting (air), layer and backing (glass) media. These refractive indices will be complex.
It is much easier to go from a situation where you know the refractive indices of everything to calculate the expected reflectivity/transmission.
In your case you have reflectivity as a function of angle/wavelength and want to calculate the refractive indices + thicknesses of the layer. Whilst you can invert the data it's fraught with problems. It's much easier to treat the problem as a curve-fitting exercise. Thus, you make a model system, with the fitting parameters being the refractive indices and thicknesses, etc. Then you calculate the expected reflectivity. Then you refine those parameters until the theoretical reflectivity matches the measured data. Differential equations or root finding is not the way to go here, trust me.
Looking at your experiment I see that you may have already coded up most of what is needed.
April 5, 2011 at 05:02 pm - Permalink
I had tried to do something along these lines but I think I got a little lost along the way. The function duplicater() in my experiment file is an attempt at finding a graphical solution (interecept of 2 curves with the xy (representing n, k in this case) gives a point, interpreted as the root). This process however encounters some strange problems, as can be seen in graphs 6 and 4. Graph 6 is a contour plot of the abs(reflectance - measured) + abs(transmission - measured) for point 49 (wavelength = 1353 nm). I add the absolute values so the minimum point is the point closest to the root. When calculating only this point, things seem to work fine. However, when I calculate points 45 through 49(updating controlled by an if-else statement at the end of the function), the program seems to be led astray by something. As seen in graph 4, the root of point 49 has shifted. The original values for n, k from graph 6 are included in the calculation which produces graph 4. As such, I am at a complete loss as to why it doesn't calculate the same root. Note: calculating more values causes the roots to move further and further away from their expected values.
As seen in the code, the calculation only uses the previous point for determining the range of n, k values to calculate. As such, the root found should be the same, as long as it exists in the range of attempted values.
Should I instead try to do a curve fitting routine (2D) of the equations to the transmission and reflectance curves over this possible range of n, k?
April 6, 2011 at 06:12 am - Permalink
Thats what I would be doing.
April 6, 2011 at 02:24 pm - Permalink
April 7, 2011 at 03:14 am - Permalink
This is something I was concerned about. It is my understanding that a curve fitting routine will give single values for n, k. Since they are wavelength dependent, am I correct in assuming that I will have to do this curve fitting routine for each wavelength?
As mentioned previously, it seems that there is no root at point 356, see the image posted by John Weeks. Will a curve fitting routine even find an answer here?
April 7, 2011 at 11:03 am - Permalink
I would like to know if you finally solved your mathematical problem. I just saw this inquary and I was wondered because I used Newton-Raphson method to solve a non-linear system equation 2X2, very similar you used. In fact, I also had to determine the complex refractive index, that is, calculate the n and k values to get m=n + ik. I you have any comment or questions feel free for asking me, ok?
Cheers
February 14, 2013 at 12:35 pm - Permalink
February 15, 2013 at 06:29 am - Permalink