Refraction equation, quartic equation

Question

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:

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<x<d$):

$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<d$. Do you have any idea how to simplify this equation or maybe you can suggest me an easier way to find $x$?

Since the solution for n=1 appears to be simple, perhaps you could let $(n^2-1)=\epsilon$ and use singular perturbations to get an asymptotic form of the solution — Valentin, May 28 '12 at 22:14
@Valentin: I need a solution for $n>1$ (or at least for $n=1.33$). I'd prefer to find an exact solution. — miloszmaki, May 28 '12 at 22:30
asymptotc expansion will give you a solution valid for some range of $n>1$ which should be computationally more efficient than horror like this. what is p? — Valentin, May 28 '12 at 22:48
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. — miloszmaki, May 29 '12 at 16:30

score 2 · Answer 1 · answered May 29 '12 at 06:57

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)$.

answered May 29 '12 at 06:57

joriki

238,052

Thank you! I will check this method if it's accurate enough. – miloszmaki May 29 '12 at 09:06
Although it generally works well for $x=dh/(h+w)$, there is a problem for this. 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. – miloszmaki May 29 '12 at 11:09
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? – miloszmaki May 29 '12 at 13:38
@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 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. – joriki May 29 '12 at 14:22
1

That's a good idea. Finding a root for $0<x<d$ is guaranteed, since $f(0)=-d^2h^2<0$ and $f(d)=d^2n^2w^2>0$ and $f$ is continous. – miloszmaki May 29 '12 at 16:24

score 1 · Accepted Answer · answered May 29 '12 at 18:12

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.

answered May 29 '12 at 18:12

For example, here's the graph for your example in the comments. I brought $n$ to the other side of the equation so you can see the solution as the zero crossing. – May 29 '12 at 18:24
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. – miloszmaki May 30 '12 at 08:47
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. – miloszmaki May 31 '12 at 19:23

Refraction equation, quartic equation

2 Answers2

Linked