はじめに
2つのベジェ曲線の交点を、陰関数を経由して代数方程式に変換し解（交点のパラメータ）を求める方法（Implicitization）について、前回、Sage Math上で解を求めました。
代数方程式の解を求める（求根）処理の実装を目指すため、DKA法を試してみました。
結果、うまくいく場合といかない場合がありました。
DKA法について
DKA法の解説や実装例は以下を参考にしました。
DKA法はニュートン・ラフソン法をベースに、高次方程式に対応した求根アルゴリズムです。複素根に対応しています。関数の次数と各項の係数が入力となります。
根を包含する円周上に等間隔で初期値を配置し、近似計算を反復します。
複数の近似解がひとつの根に収束しないよう、近似解どうしの差が漸化式の計算に反映されるしくみ（斥力の作用）が備わっています。
初期値の円の半径を求める処理が面倒くさいですが、正しく設定する必要があります。小さすぎると不安定な挙動が発生し、大きすぎると収束が遅くなります。
例①（うまくいく）
曲線0（オレンジ色）を $(x, y)=(f_0(t), g_0(t))$、曲線1（青色）を $(x, y)=(f_1(u), g_1(u))$ とします。
曲線0の陰伏方程式 $F_0(x,y)=0$ について、$x$ に $f_1(u)$、$y$ に $g_1(u)$ を代入した $u$ についての9次方程式を $F_0(u)=0$ とします。
曲線1の陰伏方程式 $F_1(x,y)=0$ について、$x$ に $f_0(t)$、$y$ に $g_0(t)$ を代入した $t$ についての9次方程式を $F_1(t)=0$ とします。
$F_0(u)=0$ と $F_1(t)=0$ の係数と根をSage Mathで求めました。それぞれ以下の通りです。
（項0は $u^9$、$t^9$ の係数、項9は $u^0$、$t^0$の係数）
| 項 | $F_0(u)$の係数 | $F_1(t)$の係数 |
|---|---|---|
| 0 | -396417457736 | 396417457736 |
| 1 | 1687229544456 | -1837788312096 |
| 2 | -2959049224008 | 3582099890784 |
| 3 | 2763498928752 | -3814300664676 |
| 4 | -1484376700032 | 2413858345488 |
| 5 | 462506013504 | -925767453096 |
| 6 | -80122864392 | 210592186686 |
| 7 | 6974963136 | -26777484144 |
| 8 | -258820272 | 1716328440 |
| 9 | 2909250 | -42600087 |
| $F_0(u)$の根 | $F_1(t)$の根 |
|---|---|
| 0.0192470307676812 | 0.0705390793450691 |
| 0.0549916382769417 | 0.0972454495384701 |
| 0.122057002218322 | 0.151317532883142 |
| 0.320989106683592 | 0.423317135673883 |
| 0.442530534075141 | 0.504784914630604 |
| 0.488155778702660 | 0.616234378772213 |
| 0.900475882753311 | 0.858627731561791 |
| 0.944674407179144 | 0.947328062423360 |
| 0.963072466334768 | 0.966598091975171 |
この係数をDKA法に入力しました。この例ではすべての根（$t$ と $u$）が $[0,1]$ の範囲の実数なので初期値の円は小さな半径で済み、迅速に収束しました。
反復の各回における近似解は以下の通りです。（左：$u$、右：$t$ ）
青点線は初期値が配置される円（半径 0.8911451625812689 と 0.8150730626130066）です。
反復の各回における近似残差（Sage Mathで求めた根との距離）は以下の通りです。
どの程度の精度が必要かによりますが、15回程度の反復で収束しました。
例②（うまくいく）
係数は以下の通りです。
| 項 | $F_0(u)$の係数 | $F_1(t)$の係数 |
|---|---|---|
| 0 | -1 | 1 |
| 1 | 18 | 27 |
| 2 | -18 | 189 |
| 3 | -6489 | -264 |
| 4 | 68040 | -3780 |
| 5 | -209412 | 7803 |
| 6 | 128169 | -2283 |
| 7 | -20250 | -5481 |
| 8 | 276534 | 10476 |
| 9 | -153711 | -3502 |
この例では、ベジェ曲線の交点は1つだけですが、実数根が5個ずつあります。
| $F_0(u)$の根 | $F_1(t)$の根 |
|---|---|
| -17.4455311309653 | -11.7491222197159 |
| 0.559001039805706 | -1.13495840453138 |
| 1.58591522923000 | 0.427184956901181 |
| 4.61782798597658 | 1.82778197153885 |
| 5.80036113872243 | 2.78366617296693 |
| -0.584112384659363 - 0.893083467241210*I | -10.1061675932554 - 1.15206211156715*I |
| -0.584112384659363 + 0.893083467241210*I | -10.1061675932554 + 1.15206211156715*I |
| 12.0253252532746 - 13.4619234074658*I | 0.528891354675582 - 0.942411682190976*I |
| 12.0253252532746 + 13.4619234074658*I | 0.528891354675582 + 0.942411682190976*I |
ベジェ曲線のパラメータの値域は $[0,1]$ ですが、それより外側を作図すると以下のようなおどろおどろしい曲線となります。交点が4つあります。自己交差もあります。
さらに外側を作図すると、はるか彼方にもう1つ交点があります。
パラメータの値域を無視すると交点は5個で、実根の数と一致します。
反復の各回における近似解は以下の通りです。（左：$u$、右：$t$ ）
初期値が配置される円の半径は 20.38309113011681 と 14.12514808245628 です。
ベジェ曲線としては、1つ以外の交点は『どうでも良い』のですが、DKA法ではそうもいかず、すべての根を包含するよう初期値の円が大きくなってしまいます。
反復の各回における近似残差は以下の通りです。
反復回数は20回程度に増やす必要があります。
例③（うまくいかない）
例②とほぼ同じですが、$P_0$ が左にずれています。
この作用で『どうでも良い』交点がはるか彼方に移動し、初期値の配置される円の半径が 624193.8162941233 と 518608.3232018586 と、とてつもなく大きくなってしまいました。
DKA法を反復50回までは試してみましたが、収束しませんでした。
ちなみに、Sage Mathで求めた係数と根は以下の通りです。
| 項 | $F_0(u)$の係数 | $F_1(t)$の係数 |
|---|---|---|
| 0 | -2.84217094304040e-14 | 2.13162820728030e-14 |
| 1 | 4.54747350886464e-13 | 1.92096649698215e-8 |
| 2 | 4.32191882282495e-8 | 0.00997905112308217 |
| 3 | -5184.01247445624 | 1727.93513568025 |
| 4 | 62208.0991688475 | -9071.84469953083 |
| 5 | -210924.160909874 | 13283.8409573164 |
| 6 | 131220.153425417 | -2375.98128942742 |
| 7 | 29888.0681867107 | -10286.8914761232 |
| 8 | 216514.307900543 | 14498.9058213240 |
| 9 | -135594.460289519 | -4589.97442830960 |
| $F_0(u)$の根 | $F_1(t)$の根 |
|---|---|
| -567115.622411213 | -248197.018538645 |
| 0.557665338289116 | -1.03330670884847 |
| 1.61807431932095 | 0.446194297282387 |
| 5.16799108368957 | 1.86453498958389 |
| 5.79814264427882 | 2.90335369440187 |
| -0.570941565478316 - 0.800872586604778*I | -213844.115726888 - 446351.819141408*I |
| -0.570941565478316 + 0.800872586604778*I | -213844.115726888 + 446351.819141408*I |
| 283559.811210479 - 491136.142417216*I | 0.534608073955102 - 0.882297985704637*I |
| 283559.811210479 + 491136.142417216*I | 0.534608073955102 + 0.882297985704637*I |
反復回数を大幅に増やせば収束するかもしれませんがそもそも、9次関数その他の計算誤差がひどいことになっていると思われます。64ビットdoubleでは厳しいかと。
また、これがワーストケースというわけではありません。