4
$\begingroup$

Given two points $P$ and $Q$, a line ($A$, $B$ - orthogonal projection of $P$, $Q$ onto the line) and a coefficient $n$, I want to find out such point $C$ that $\frac{\sin{a}}{\sin{b}}=n$ (in fact, it's an equation of light refraction). I also assume that $C$ lies between $A$ and $B$. See the image below:

image

What I know is: $w, h, d, n$. What I want to calculate is $x$.

I need to find the easiest solution possible, in order to calculate $x$ in a computer's application very efficiently.

What I've found out so far is:

$$n=\frac{\sin{a}}{\sin{b}}=\frac{\frac{d-x}{\sqrt{(d-x)^2+w^2}}}{\frac{x}{\sqrt{x^2+h^2}}}=\frac{(d-x)\sqrt{x^2+h^2}}{x\sqrt{(d-x)^2+w^2}}$$

Hence, we need to solve this equation for $x$ (assuming $0

$n^2x^2((d-x)^2+w^2)=(d-x)^2(x^2+h^2)$

which yields:

$f(x)=(n^2-1)x^4-2d(n^2-1)x^3+((d^2+w^2)(n^2-1)+w^2-h^2)x^2+2dh^2x-d^2h^2 = 0$

That way we obtained a quartic equation. How to solve it? I've tried using the Ferrari's solution, but the result is very complicated. Moreover, this equation can have as many as four real roots but I know it's got only one root for $0$x$?

  • 0
    Since the solution for n=1 appears to be simple, perhaps you could let $(n^2-1)=\epsilon$ and use [singular perturbations](http://en.wikipedia.org/wiki/Perturbation_theory#Example_of_second-order_singular_perturbation_theory) to get an asymptotic form of the solution2012-05-28
  • 0
    @Valentin: I need a solution for $n>1$ (or at least for $n=1.33$). I'd prefer to find an exact solution.2012-05-28
  • 0
    asymptotc expansion will give you a solution valid for some range of $n>1$ which should be computationally more efficient than horror like [this](http://www.wolframalpha.com/input/?i=%28n2−1%29x4−2d%28n2−1%29x3%2B%28%28d2%2Bw%5E2%29%28n2−1%29%2Bw%5E2−p2%29x2%2B2d*h%5E2x−%28d%5E2%29*h%5E2%3D0). what is p?2012-05-28
  • 0
    $p$ was supposed to be $h$; I fixed it.2012-05-29
  • 0
    @joriki: Right, thanks for fixing it ;)2012-05-29
  • 0
    Have you heard of Snell's Law?2012-05-29
  • 0
    Yes, but deriving $x$ from Snell's Law gives the same result - quartic equation. In fact, my problem is to find the point of refraction of light according to Snell's Law.2012-05-29

2 Answers 2

1

As long as you're using a numerical method, you might as well stick with a non-polynomial form of your equation that is more well-behaved. In particular, you could just solve $$\frac{(d-x)/\sqrt{(d-x)^2+w^2}}{x/\sqrt{x^2+h^2}}=n$$ for $x$. The left-hand side is a monotonically decreasing function of $x$ between $0$ and $d$, because the numerator is nonnegative and decreasing and the denominator is nonnegative and increasing. So you can just apply bisection search and it is guaranteed to work.

  • 0
    For example, [here's the graph for your example in the comments](http://www.wolframalpha.com/input/?i=plot%20%28%28d-x%29/sqrt%28%28d-x%29%5E2%2bw%5E2%29%29/%28x/sqrt%28x%5E2%2bh%5E2%29%29%20-%201.33%20where%20w=1,%20h=5,%20d=30%20for%20x=0..30). I brought $n$ to the other side of the equation so you can see the solution as the zero crossing.2012-05-29
  • 0
    Right, converting the original equation into a polynomial introduces unwanted ambiguity. The polynomial can have as many as four roots while there's only one correct solution to the original problem. I just hoped that I could find a descent analytic solution but it turns out that a numerical method is the only way to solve the problem efficiently.2012-05-30
  • 0
    I've tested it and it works great! Only 10-20 iterations of bisection search is enough. I could also improve the solution using Newton's method.2012-05-31
2

A plot of a typical example shows that the right-hand side of your last displayed equation behaves rather nicely and it should be possible to obtain the solution efficiently using Newton's method, which wouldn't require drawing any roots and should yield the result to sufficient accuracy within a few iterations that only involve a couple of multiplications and additions and one division each.

You could use the solution $x=dh/(h+w)$ for $n=1$ as the initial value for the iteration.

By the way, I'd divide through by $d^4$ to get rid of the irrelevant scale and express everything relative to $d$:

$$\epsilon\xi^4-2\epsilon\xi^3+((1+\omega^2)\epsilon+\omega^2-\eta^2)\xi^2+2\eta^2\xi-\eta^2=0\;,$$

with $\epsilon=n^2-1$, $\xi=x/d$, $\omega=w/d$ and $\eta=h/d$ with initial value $\xi_0=\eta/(\eta+\omega)$.

  • 0
    Thank you! I will check this method if it's accurate enough.2012-05-29
  • 0
    Although it generally works well for $x=dh/(h+w)$, there is a problem for [this](http://www.wolframalpha.com/input/?i=plot+%28n%5E2-1%29+x%5E4+-+2+d+%28n%5E2+-+1%29+x%5E3+%2B+%28%28d%5E2%2Bw%5E2%29%28n%5E2+-+1%29+%2B+w%5E2+-+h%5E2%29+x%5E2%2B2+d+h%5E2+x+-+d%5E2h%5E2+where+n%3D1.33%2Cw%3D1%2C+h%3D5%2Cd%3D30). Starting with $x=25$ the result is getting close to $x=30$ while $x=30$ is not a root. The correct result can be obtained starting with $x=0$ which gives approx. $5.7$. Also I'm not sure if starting with $x=0$ always gives the correct result.2012-05-29
  • 0
    Edit to my previous comment: when starting with $x=0$, it converges to $-5.7$ (which means $n=-1.33$). The correct result can be obtained starting for example with $x=1$ or $x=12.5$. Do you have any idea how to choose the starting point which guarantees the convergence and the correct result?2012-05-29
  • 0
    @miloszmaki: No, sorry, I wasn't expecting such difficult behaviour in other cases from the example I tried. You could use more robust methods like [regula falsi](http://en.wikipedia.org/wiki/Regula_falsi) instead (starting with $x\in[0,d]$), or use such a method to get close and then switch to Newton's method to improve convergence.2012-05-29
  • 1
    That's a good idea. Finding a root for $0 is guaranteed, since $f(0)=-d^2h^2<0$ and $f(d)=d^2n^2w^2>0$ and $f$ is continous.2012-05-29