Rearranging a Function to Find X

Thanks mate, nice one. Is that MathCAD? - I don't have .

Actual values are a bit different to your assumptions:

0.5 < x<5 (so no need to fret about the singularity),

A ≈ 10

B ≈ 0.001

The exponential term is just a small correction (to what's basically an inverse square law), but it's important to the end user.

I tried taking logs, but (as you found) it gets just as messy. Got a horrible feeling I'll end up with some kind of interative solution, after setting B to zero to get a first approx.

Yes, it's Mathcad.

Is there a particular value for D you are trying to get? You imply that you know the value of D(x), and are trying to find the value of x for that D(x). Correct?

I know the value of D for one specific value of x. The object is to be able to calculate x for any required value of D.

Might as well fill you in a bit more.

This is for use on a rig for calibrating radiation detectors. One of a number of radiation sources is exposed, producing a gamma ray beam which shines onto the detector which is mounted on a motorized trolley. Once a year or so, the doserate in the beam is measured at one specific point (the "calibration distance") using a secondary standard ionization chamber. This gives what's known as the "source term".

To be able to expose the detector to a specific doserate, the source term is first corrected for decay (using a simple calculation based on the half-life of the nuclide), then the required distance from the source to the dectector must be calculated.

The doserate falls off in the usual (inverse square law) manner, but there is also a small amount of attenuation due to photoelectric and Compton interactions between the gamma photons and the air.

The full expression for doserate is:

D(x) = Dc * (Xc/x)² * e^-μ(x-Xc),

where:

Dc is the decay-corrected "source term",
Xc is the calibration distance, and
μ is the air attenuation coefficient.

An interesting quirk is that for some low energy radiations, the sign of μ can change, making the resulting doserate higher than would be expected from just the inverse square law calculation.

In use, an operator selects a nuclide and types in a required doserate (or these are found from a protocol file for a particular instrument to be calibrated). The software then finds the most appropriate combination of radiation source and distance, moves the trolley to the correct position, moves the source into position and opens a shutter to start the exposure.

Many thanks for your help - including the pm's that others haven't been privy to.

_{(Jeez - I've ended a sentence with a preposition again! I'll be in trouble with the Grammar Police!)}

By the way, if B = 0.001 then across the range from x = 0.5 to x = 5, the e^-Bx term only changes value from 0.9995001 to 0.9950125.

So you could get an excellent approximation by simply setting the e^-Bx term to K = 0.997 [or even just to K = 1.00]. Then solve D(x) = Ax^-2·K for x. Easy as pi.

Agreed - thing is that the B can have other values, and does just about become significant in some circumstances.

I was considering using a linear approximation (which would give a result much more accurate than many of the other uncertainties in the system), but thought it would be more elegant if I could show that the software's doing a "proper" calculation.

MATHCAD(R), like other programs, can be powerful analytic tools, but maybe Excel(R) can help out until someone provides a more exact formula to answer your question. In the mean time, even though I didn't go to a theory-heavy school, maybe this can tide you over 'til then.

The basic equation is, essentially,: d = A(x^(-2)) * e^(-Bx), such that:

A = 10; B = .001; and 0.5 < x < 5.0 => one dependant, & one independant, variable

=> d solutions should be unique for each variable x => take natural logs of both sides

=> ln(d) = ln (Ax^(-2)) + ln( e^(-Bx) = ln(A) - 2 ln(x) - Bx

=> which yields quantifiable equations of the form;

ln(d) = ln(A= 10) - 2 ln(x = 0.5) - (B = .001) (x =0.5)

.

.

.

ln(d) = ln(A= 10) - 2 ln(x = 5.0) - (B = .001) (x = 5.0),

from which unique numerical values for e^(ln(d)) = d(x), can be calculated, which readily allow for interpolation to get the intermediate x-values.

This furthur suggests using Excel(R) to plot-out d(x)-vs-x, or better yet, to plot x-vs-d(x), which provides a graph from which to read-off the x-value from the d(x) curve, as was asked for in the original posting ...

Don't really understand why you bothered to take logs - the exponential function is available in Excel and also in most (if not all) programming languages.

The following isn't directly in reply to you, Guest (but thanks for the input anyway) - I'm just thinking.

Your post has suggested to me a very simple method of finding my answer.

• I know that the smallest increment I need for x is 0.0001 (a consequence of the electromechanical system controlling the distance x from the origin).
• Taking steps of x of 0.0001 in the range 0.5 < x < 5.0 gives 45,000 values.
• All 45,000 values for D(x) could be computed in well under a second on a PC.
• The function D(x) is well-behaved in this range.
  
• Given this array of values of D(x), it would be very easy to find the appropriate value of x for any value of D between D(0.5) and D(5.0).

This may be the solution I adopt!

My muse is this:

Is the accessibility to "number-crunching machines" leading to a neglect of the study of mathematical methods?

We each approach a task with the particular tools which we have at hand at the time, and Excel(R) is what I have...

In any case, I am glad that I was able to be of assistance. Happy New Year, in any case...

I hate to be critical, but here you would use ln, not log.

Still, GA for the input and the workthrough on Mathcad. I'd like to get a license, but can't really justify the expenditure right now.

True, but really, 'log' is generic. Though ordinarily 'log' is considered to be 'log₁₀' and 'ln' is 'log_e', I didn't make the distinction since the equation had the e in it. I can see how someone might be confused since I was a bit causal about that. So, good point.

I also tend to be a bit sloppy with the use of the word. I'm writing the code in VB, and the log() function actually returns the natural log.

_{(That's my excuse anyway ).}

Just a quicky... because we are working with e^x be sure to use ln (natural log) rather than log₁₀.

Bill

Given that the unknown is both inside and outside a log, there is no expression directly for x.

Since it is being done on a PC, numerical methods are your best option! I don't remember them all, but I do remember that the Runge-Kutta (sp?) method was a pain to program, but is extremely efficient.

Again, no direct solution, but I just remembered a quick and simple numerical method... for your equation it would be the following

x_n+1 = (1/B)*ln[D/(y*x²)]

Just plug in your particular D and a guess for x on the right side to compute the result. Use the result to plug back into the right side again... repeat. It will converge relatively quickly depending on the desired level of accuracy.

Hope it helps!

Thanks - I'll give it a try (if I can work out what 'y' is ).

My apologies, I was in a rush on my last post. Here's the correction and the method.

x_n+1 = 1/B*ln[ A / (D*x²) ]

If you take your original function and simply rearrange to solve for one "x" you can obtain an iterative solution to the function. That's all the above amounts to. To use it you simply enter in your known D(x) value and a guess for x on the right, compute the answer, and use that solution for the next iteration.

It is relatively quick (about half of the error is reduced each iteration).

Thanks again - will try it.

Just to save me having to think (still a bit foggy after sampling some of the things Santa's brought!) - I assume that's

x_n+1 = [1/B]*ln[ A / (D*x²) ] ?

I was considering using a binary chop, or even a linear search - I know the attenuation corrected solution's within a couple of mm of the inverse square result, and I only need a result to 0.1mm, but this looks nice.

You have it correct.

There very well may be other solutions/methods that would converge faster. This was a "quick and dirty" method a professor taught me. The method doesn't always work, but in this case I tried it out in Excel and seemed to work great.

Earlier it sounded like the OP was worried about having an iterative solution. He also would like the solution to be elegant. Since some people seem to be unaware of iterative solutions, an iterative solution IS elegant.

I did use the comparative - more elegant.

I hope this doesn't appear too silly and please correct me if I'm wrong but wouldn't the rules of simple algebra indicate that dividing both sides of the equation by D is in order here?

D(x) = Ax^-2.e^-Bx

(D(x))/D = (Ax^-2.e^-Bx)/D

x = (Ax^-2.e^-Bx )/D

This is a differential equation. I'm nowhere near my books, and it's been a little while, or I could give you a better answer, but there is no real simple solution to your problem. I can't remember which answer it is, but here are your two options:

1) Solutions of the form x = Asin(y) + Bcos(y)

2) This may be a transcendental equation. What that means is you can't actually solve it.

A linear approximation may be your shortest route to victory.

try this ones

y = A/x^2 - AB/x + 1/2 A B^2 -1/6 A B^3 x + 1/24 A B^4 x^2 - 1/120 A B^5 x^3

this is an approximation of the orginal function.

e^x = Σ xⁿ /n!

e^x = 1 + x/1! + x²/2! + x³/3! + x⁴/4! + x⁵/5! + ….

e^-Bx = 1 - Bx + B²x²/2 - B³x³/6 + B⁴x⁴/24 - B⁵x⁵/120 + ….

Ax^-2 e^-Bx = Ax^-2 - ABx^-1 + AB²/2 - AB³x/6 + AB⁴x²/24 - AB⁵x³/120 + ……

If you take A = 10 and B = 0.001 and x = 0.5 then increases the function with about minus 0.00000000083 if you take only the first five terms

In the region of x = 0.5, the increase of the function in maximal when x changes.

If this is still not accurate enough, you have to take more terms.

Thanks, Rudy - now solve for x .

Seriously, your expansion is useful for plugging into the Newton-Raphson method for finding the solution by iteration.

Just taking the first term of your expansion, with f(x) = (Ax^-2 e^-Bx) - D,

f(x) ≈ Ax^-2 - D,

f'(x) ≈ -2Ax^-3

which is what I used in my solution in #28.

Including the second term would give me

f(x) ≈ Ax^-2 -AB^x-1 - D, and

f'(x) ≈ -2Ax^-3 + ABx^-2

I'll try it for completeness (?) - but I don't think it will make a practical difference. Thanks anyway.

Hey john, thanks for the GA

What is the error, the fault tolerance you can accept ?

Which are the smallest values of D(x) and x you want to calculate ?

For x it is about 0.0001 according your post number 27, but what about D(x) ?

How accurate are the values of A and B ?

Rudy

Hi, Rudy,

As I've said in previous posts, the full equation is:

D(x) = Dc * (Xc/x)² * e^-μ(x-Xc),

where:

Dc is the decay-corrected "source term",
Xc is the calibration distance, and
μ is the air attenuation coefficient.

0.5000 <= x <= 5.0000, with x determined to the nearest 0.0001

Xc = 1.0000 (usually, though it may sometimes be 1.5000 or possibly 2.0000).

Dc can be anything from 1.000E0 to 1.000E-6 (it depends on the radiation source).

The range of D(x) is limited by the range of x values; with Xc = 1 and Dc normalized to 1 the range would be approximately D(0.5) = 4 to D(5) = 0.04. Note that Dc is just an abreviation of D(Xc).

Dc and μ are accurate to 4 significant figures.

The error in the answer should be as small as possible, but it is limited by the fact that x can only take the discrete values as noted above.

I've found something that works!

From my previous post #10, the full expression for doserate is:

D(x) = Dc * (Xc/x)² * e^-μ(x-Xc),

where:

Dc is the decay-corrected "source term",
Xc is the calibration distance, and
μ is the air attenuation coefficient.

On the basis that the exponential part changes very slowly compared with the power part, I thought I'd try a slightly "bent" version of the Newton-Raphson method.

Let f(x) = Dc*Xc²*x^-2 * e^-μ(x-Xc) - Dr,

where Dr is the required doserate.

A good approximation for the first derivative is then

f'(x) = -2Dc*Xc²*x^-3

Plug this into the Newton-Raphson expression:

x_n+1 = x_n - f(x_n)/f'(x_n)

Using the inverse square law result as a first approximation for x, this converges to give a value of x yielding the required doserate within 0.01% after 4 iterations, even with |μ| as large as 0.1.

Many thanks to all for your invaluable input, which got me thinking hard enough and in the right direction to get there.

I'm glad you found a solution. I love these types of problems and had a great time playing with it. I tried out Newton's method as well and probably should have communicated my findings (-1 point for me). I found Newton's slower or found a different solution outside the range you're interested in than the "quick and dirty" method (and hence I should have inquired more).

But that aside, good job and thanks for the brain exercise!

I tried your solution (in #22) using Excel, and my attempt diverged violently! Think I must've misinterpreted something somewhere along the line. I'd be interested to see your method working - if you're OK with sending me an xls, please pm me and I'll give you my e-mail address. I'll understand if you're reluctant.

Thanks anyway for the prod in the right direction. I've found that for realistic values of the attenuation coefficient, my bodge converges sufficiently within 2 iterations, so I'm pretty happy with it. Provided I start with the (easily calculated) inverse square result as a first approximation, I haven't managed to find any "wrong" solutions.

Glad you had fun!