View Question
Q: example of iterative method to solve equation ( Answered ,   0 Comments )
 Question
 Subject: example of iterative method to solve equation Category: Computers > Algorithms Asked by: slytech-ga List Price: \$200.00 Posted: 21 Apr 2006 10:10 PDT Expires: 21 May 2006 10:10 PDT Question ID: 721381
 ```I am trying to implement a program that uses Allard's Law to calculate the visual range (R) of lights at night. The formula I have found is T = (I*e^(-1*x*R))/(R^2), where T is the visual threshold of illumination, x is the extinction coefficient, and I is the intensity of the lights. From what I've read, the formula cannot be solved for R directly, and one has to use "iterative methods" to to determine it. This is where my knowledge of mathematics break down. Question: Demonstrate using a plain C/perl/basic/pseudo-code example, how I can best determine R given the value for the other variables. It needs to be as complete, simple and accurate as possible.```
 ```Dear slytech, When a formula does not have a closed-form solution, an iterative approach such as Newton's Method can often be used to arrive at an arbitrarily precise answer by computing a sequence of successively closer approximations. Arbitrary precision means that the user can set a margin of error as small as he likes and then carry out the iteration until it produces an answer within the margin of error. In some cases, iterative approximation doesn't work because the formula converges toward false local solutions or diverges altogether in some regions of the curve. Fortunately, Allard's Law is proven to converge toward the true solution for all values that are meaningful in the physical world. The Aerodrome Meteorological Observing Systems Study Group (AMOSSG) of the International Civil Aviation Organization describes in its Manual an efficient way to iteratively solve Allard's Law for the variable you call R. See page 17 of the following PDF document. Draft AMOSSG Manual: Appendix B [pdf] http://www.icao.int/anb/SG/AMOSSG/meetings/amossg4/sn5B.pdf Note that AMOSSG uses different variable names than yours, but otherwise gives exactly the same formula for Allard's Law. What you call T, the visual threshold of illumination, is E_T in the Manual; what you call x, the extinction coefficient, is sigma (Greek lowercase s) in the Manual; and R, the visual range or visibility, is called V in the manual. The name of variable I, for luminous intensity, is unchanged. Although AMOSSG mentions that I should be 1000 cd when calculating visibility for aerodrome purposes, you are of course free to use any other value without affecting the validity of the formula or its iterative solution. Also note that the variable MOR, which stands for Meteorological Optical Range, is defined on page 16 as 3/s, so it really doesn't introduce anything new. I shall quote the relevant passages from page 17 before giving you my code example. I have made some minor typographical corrections to the text. Call E_T the visual threshold of illumination and I, the luminous intensity. With E_T = I * e^(-s*V) / V^2 and by replacing s by 3/MOR, we get: V = -MOR/3 * ln(E_T / I * V^2) (1) [...] Consider the sequence V_n = -MOR/3 * ln(E_T / I * V_{n-1}^2) = f(V_{n-1}) . If this sequence converges, it converges toward V, the visibility sought. It can be demonstrated that if V_n is greater than V, then V_{n+1} will be less than V. The sequence V_n is close to the value of V. If we take V_0 = MOR and if V_1 is less than V_0, we can conclude that solution V to equation (1) is less than V_0 = MOR. In this case, the calculation is not necessary, since the visibility distance given by the law of Allard is less than the MOR. So the visibility is equal to the MOR, which is good since the sequence can diverge in such conditions. However, it is possible to show that if V_1 > MOR, given that V_0 = MOR, the sequence converges. The iterative calculation of this sequence can then be performed until the difference between V_n and V_{n+1} is small in relation to the value of V_n. For example: abs(V_n - V_{n-1}) / V_n < 0.01 . In practice, convergence may be slow. It can be very quickly accelerated using an intervening variable: Start with V_0 = MOR and calculate V_1 = f(V_0). Calculate V_01 = (V_0 + 2*V_1) / 3. Calculate V_2 = f(V_01) and V_12 = (V_01 + 2*V_2) / 3. Calculate V_3 = f(V_12) and V_23 = (V_12 + 2*V_3) / 3. The calculation can be continued so forth, but in practice, the value of V_23 is very close to the required value of V and the calculation can be ended at the third iteration. In the C code below, I use your variable names for T and x, but the Manual's names for elements in the V_n sequence. We are solving for what you call R, which the Manual calls V. double f(double V, double T, double I, double MOR) { return -MOR/3 * log(T / I * V*V); } double compute_visibility(double T, double x, double I) { double MOR, V_0, V_1, V_01; if (x == 0) { printf("error: extinction coefficient must not be zero\n"); return -1.0; } V_0 = MOR = 3.0/x; V_1 = f(V_0, T, I, MOR); if (V_1 < V_0) return MOR; V_01 = (V_0 + 2*V_1) / 3.0; while (abs(V_1 - V_0) / V_1 >= 0.0001) { V_0 = V_1; V_1 = f(V_01, T, I, MOR); V_01 = (V_01 + 2*V_1) / 3.0; } return V_1; } The margin of error is currently set to 0.0001, or 0.01%. You may adjust this value in the test condition of the while loop. I have enjoyed answering your question. If you have any concerns about the completeness or accuracy of my research, please advise me through a Clarification Request and give me a chance to fully meet your needs before you rate this answer. Regards, leapinglizard Search strategy: allard law range illumination ://www.google.com/search?hs=gWJ&hl=en&lr=&client=firefox-a&rls=org.mozilla%3Aen-US%3Aofficial&q=allard+law+range+illumination&btnG=Search allard law iterative solution ://www.google.com/search?hs=zqy&hl=en&lr=&client=firefox-a&rls=org.mozilla%3Aen-US%3Aofficial&q=allard+law+iterative+solution&btnG=Search allard law approximation range illumination ://www.google.com/search?hs=fSV&hl=en&lr=&client=firefox-a&rls=org.mozilla%3Aen-US%3Aofficial&q=allard+law+approximation+range+illumination&btnG=Search``` Request for Answer Clarification by slytech-ga on 24 Apr 2006 02:50 PDT ```leapinglizard, thanks for your quick response and sample program. Some questions: 1. I don't quite understand what "convergence" means wrt to the formula. It would be great if you could provide a few good links to help me understand the solution better. 2. What is the significance of the "if (V_1 < V_O) return MOR;" portion of the example? Is it related to the scenario the AMOSS example gives, i.e. para 4.3.6 mentions "... In this case, the calculation is not necessary, since the visibility distance given by the law of allard is less than the MOR". Or is it related to the solution, in which case the example can only be used when one is not interested in V when visibility is greater than MOR. 3. I've tried to run the example provided and compare it with some other sample data. The results match exactly, however it seems to break whenever "V_1 < V_O" is true (plus a few other instances). If it was necessary to calculate V under such conditions, can the example be improved, or is there an alternative solution? I've posted the sample data and my results using the example program you provided at this URL: http://ja.slytech.googlepages.com/home. The sample data shows the results for V for a number of scenarios, e.g. when MOR is 10000, 3000, 1000, 300, 100 and 30, for E_t of 10^-4 and 10^-6, and and I of 10000, 1000, 100 and 10. Thanks!``` Clarification of Answer by leapinglizard-ga on 27 Apr 2006 14:03 PDT ```1. I don't quite understand what "convergence" means wrt to the formula. It would be great if you could provide a few good links to help me understand the solution better. Convergence just means that the sequence gets closer and closer to a single value as you go further along. Some sequences converge for all values; some converge only for certain values, which is called conditional convergence. The question of whether a sequence of a given form converges, and under what conditions, is a thorny one that calls for deep mathematics. There is, in general, no simple way of deciding this question. Mathworld: Convergent series http://mathworld.wolfram.com/ConvergentSeries.html Mathworld: Convergence Improvement http://mathworld.wolfram.com/ConvergenceImprovement.html Mathworld: Limit http://mathworld.wolfram.com/Limit.html 2. What is the significance of the "if (V_1 < V_O) return MOR;" portion of the example? Is it related to the scenario the AMOSS example gives, i.e. para 4.3.6 mentions "... In this case, the calculation is not necessary, since the visibility distance given by the law of Allard is less than the MOR". Yes, your first hypothesis is correct. My code was a direct implementation of the AMOSS instructions. 3. I've tried to run the example provided and compare it with some other sample data. The results match exactly, however it seems to break whenever "V_1 < V_O" is true (plus a few other instances). If it was necessary to calculate V under such conditions, can the example be improved, or is there an alternative solution? You're quite right about the discrepancies. I may have misunderstood the AMOSS instructions, or perhaps they are meant for a limited set of circumstances. I find that if I swap the initial V_0 and V_1 values to ensure that they are in increasing order, my results agree with those in the table for all except the six lower-left elements. 4839 / 13400 2653 / 5722 1340 / 2468 572 / 935 247 / 373 94 / 133 -1745 / 8646 1497 / 4090 865 / 1881 409 / 749 188 / 309 75 / 113 -2921 / 4839 -474 / 2653 484 / 1340 265 / 572 134 / 247 57 / 94 -7732 / -1745 -917 / 1497 -174 / 865 150 / 409 86 / 188 41 / 75 My revised code is in the following file. http://plg.uwaterloo.ca/~mlaszlo/answers/allard.c I've been thinking of a different way to solve Allard's Law that works for all solutions within a reasonable range, such as [10^-9, 10^9]. Let me get back to you. leapinglizard``` Clarification of Answer by leapinglizard-ga on 27 Apr 2006 15:23 PDT ```I have written a new program that should work as long as the true solution is between the constant values LOW and HIGH. The program subdivides this range into a number of intervals given by the constant INTERVALS and determines how closely the midpoint of each interval matches the formula. The endpoints of the closest matching interval become the new low and high values of the entire range. This procedure is repeated a number of times given by the constant ITERATIONS. 4839 / 13399 2653 / 5721 1339 / 2467 571 / 934 246 / 372 93 / 132 2255 / 8645 1496 / 4090 864 / 1881 408 / 748 187 / 308 74 / 112 877 / 4839 703 / 2653 484 / 1339 265 / 571 133 / 246 57 / 93 302 / 2255 275 / 1496 225 / 864 149 / 408 86 / 187 40 / 74 http://plg.uwaterloo.ca/~mlaszlo/answers/allard_beta.c Using a value of 100 for INTERVALS means that with each iteration, the size of the range is multiplied by 10^-2. Doing this ten times makes for a total factor of 10^-20, which means that the original range of 10^9 is shrunk down to 10^-11. Thus, my approximated solution is within 0.00000000001 of the true solution. I believe the textbook table is off by 1 in several cases. This may be due to rounding errors or to an inadequate number of iterations on their part. leapinglizard```