How to avoid RATNZ when NSolving an equation in mathematica

Question

I Need to solve the following non-reducible transcendental equation. I have done as follows:

Data = Import["FILE NAME", "Table"]; (*get actual data from a file*)
l = Range[31, 249];
(*generate G,t pairs for the equation*)
t = Data[[l, 1]]*10^-3; (*time in ms, convert to s*)
G = Data[[l, 2]];
k = (4*Pi*1.33/633*^-9)*Sin[10 Degree/2];
sb = 0.25*^-3;(*beam size in m*)
\[Omega] = 2;(*speed of rotation in rad/s*)
s = \[Omega]/(2 Degree);
\[Sigma] = Table[0, {Length[l]}]; (*data holder, to be filled in for loop*)
For[j = 1, j <= Length[l], j++,
  y = sig /.NSolve[{G[[j]] == Exp[-((k*s*sb*t[[j]])^2)/(2*(k*sb*sig*t[[j]])^2 + 1)]*1/Sqrt[2*(k*sb*sig*t[[j]])^2 + 1]), sig >= 0}, sig, Reals][[1]];
  \[Sigma][[j]] = {t[[j]]*10^3, G[[j]], y};
  ];
\[Sigma]

Everything works until the actual attempt to find the values for sig. Then I begin getting messages:

Solve::ratnz: Solve was unable to solve the system with inexact coefficients. The answer was obtained by solving a corresponding exact system and numericizing the result.

which causes the results to become very inexact. For example at one point I accidentally put a 2 in the numerator of the exponential, and this had no effect on the output results. This means this Ratnz procedure is heavily approximating the solution, to the point of near insensitivity to the actual inputs. As you can see from the following example data, the data is not so different on a point by point basis that insensitivity to a factor of 2 is allowable.

Example Data (will be read in from text file during import):
t(ms)               G(t)
  1.00000E-003    8.09982E-001
  1.20000E-003    8.05885E-001
  1.40000E-003    8.02226E-001

How do I get mathematica to solve the system at hand, not a very poor approximation thereof?

score 0 · Answer 1 · answered Nov 21 '13 at 17:59

Since numerical precision issues are very dependent on the actual numbers used and you haven't shown two or three example lines from your data file that really highlight the problem it is then almost impossible to suggest a solution and completely impossible to actually verify that the solution will work with your data file.

Pick out two or three lines from your data file that really demonstrate the problem well, put those in an added comment and then I'll take a few minutes and see if I can find a way to minimize or perhaps even fix the problem.

You also seem to have what may be a scrape-n-paste error in your code. There may be an extra ) at the end of your first item in your list inside NSolve, either that or there is an even bigger problem that this is just hinting at.

Mathematica, unlike FORTRAN, keeps track of the precision of every value being used and then propagates that precision through all the calculations. So if you introduce a number with one digit of precision and then do calculations with that and a dozen different other values, some with more precision, you can easily end up with a result with less than one digit of precision. Since you have used constants with only two digits of precision and your calculation may need substantially more to get a good result this may be part of your problem, but we can't tell without the data.

The data points are very helpful. Thank you. Could you provide, by any means you wish, what a desired solution would be for one of those points? That way I could compare the results of my changes to your goal. — Bill, Nov 21 '13 at 23:33
If I wrap SetPrecision[Import[...],30] to force input to have 30 digits precision the result doesn't change. If I change Pi*1.33 to Pi*4/3 and 0.25*^-3 to 1/4/1000 the results change substantially, but I don't know if these changes are exactly correct. — Bill, Nov 21 '13 at 23:49
You can also look at adding WorkingPrecision->30 as an option to your NSolve. Look at the help page for NSolve, click on "Details and Options", read the info on this and then down under examples click on Options and then click on WorkingPrecision and look at that. Change any approximate decimal values to exact rationals if that is what the value really is, run your code again and see if the result is closer to the answers you know you should get. — Bill, Nov 22 '13 at 00:01
I don't know what the answers are ahead of time. I am trying to solve for them. This is a case of the answer not being known ahead of time for any data. — Elliot, Nov 22 '13 at 17:11

Bill · Answer 2 · 2013-11-22T17:15:43.523

0

Manually eliminate all decimal points to circumvent Mathematica's automatic machine precision approximation process.

In[1]:= Data={{10/10^4,809982/10^6},{12/10^4,805885/10^6},{14/10^4,802226/10^6}};
l = Range[1, 3];t = Data[[l, 1]]*10^-3;G = Data[[l, 2]];
k = 4*Pi*133/100/(633*10^-9)*Sin[10 Degree/2];
sb = 25/100*10^-3;\[Omega] = 2;s = \[Omega]/(2 Degree);
\[Sigma] = Table[0, {Length[l]}];
For[j = 1, j <= Length[l], j++, 
  y=sig/.NSolve[{G[[j]]==Exp[-((k*s*sb*t[[j]])^2)/(2*(k*sb*sig*t[[j]])^2+1)]*1/
  Sqrt[2*(k*sb*sig*t[[j]])^2+1], sig>=0}, sig, Reals, WorkingPrecision->50][[1]];
  \[Sigma][[j]] = {t[[j]]*10^3 + 0., G[[j]] + 0., y};
];\[Sigma]

Out[7]= {
{0.001, 0.809982, 888.07100020205412013305651535199097225624483722309},
{0.0012, 0.805885, 750.32195491553116360487109372989348411883522775346},
{0.0014, 0.802226, 650.84411888264632019635531363713286343974593136657}}

No ratnz warnings, no approximations (other than the 50 digits specified inside NSolve)

Compare with the results from your original code

Out[]= {
{0.001, 0.809982, 888.071},
{0.0012, 0.805885, 750.322},
{0.0014, 0.802226, 650.844}}

edited Nov 22 '13 at 17:15

answered Nov 22 '13 at 16:58

Bill

3,664
1
12
9

Not doable. This input is from an instrument, so the formatting of the data in the text files is not under my control. Manually editing it to allow for this is not practical. Also your time units are not equivalent. – Elliot Nov 22 '13 at 17:09
Sorry for the error in time units, my mistake in formatting. Corrected. Learn how to use string replacement to remove the decimal points and scale to provide exact results. I will fiddle to do that in a few minutes. – Bill Nov 22 '13 at 17:18
string = "1.00000E-003 8.09982E-001"; FullForm[StringReplace[string, {"E-003"->"*10^-6", "E-001"->"*10^-6", "."->""," " -> ","}]] can do your conversion from decimal to exact and ToExpression then turns strings into numbers. – Bill Nov 22 '13 at 17:31
You can either use SetPrecision, which appends "binary zeros" to a floating point value, or StringJoin to append "decimal zeros" to a string form and then convert to floating point, to get the higher precision needed to give you the more exact solutions. Google for Mathematica Import scientific form to see examples of how other stackoverflow questions have handled this. – Bill Nov 22 '13 at 17:56

How to avoid RATNZ when NSolving an equation in mathematica

2 Answers2