0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

陰関数でベジェ曲線の交点を求める ~ その3:ブダン・フーリエの定理とヴァンサンの定理で根分離しリダーズ法で根を求めてみた

Posted at

前回、9次方程式の根をDKA法で求めようとしましたが、ベジェ曲線のパラメータの値域(0~1)から大きく外れた根(不要な根)に邪魔されて目的とする根を求めることが困難、との課題に直面しました。

今回、NVIDIAのAlexander Reshetovさんの論文を参考に、ブダン・フーリエの定理とヴァンサンの定理で所定の値域の実根を分離し、リダーズ法で根を求める方法を試してみたところ、期待通りの結果を得ることができました。高次方程式の根を合理的に求める方法がわかったので、陰関数を用いてベジェ曲線の交点を求める手法が実装できそうです。

参考にした論文について

多数のベジェ曲線でモデリングされたPBRTのhairやfurを対象としたレイトレーシングについて;

  • 3次ベジェ曲線と光線(直線)の交点を求める。ベジェ曲線を太さのあるものとして扱う必要があり、曲線と直線の間隔の2乗(6次関数)が最小になる位置を求めるため、その微分(5次関数)の根を求める
  • 根を求めるため、ブダン・フーリエの定理により絞り込んでヴァンサンの定理により実根を分離し、求根アルゴリズムをなるべく反復なしで適用する
  • 求根アルゴリズムはニュートン・ラフソン法、セカント法、ミュラー法、リダーズ法などを比較
  • ベジェ曲線と直線の空間的な関係性から根の存在する範囲を絞り込む、等のテクニックを組み合わせることにより、リダーズ法反復なし1回とセカント法の組み合わせで十分な近似精度が得られる(リダーズ法2回はやりすぎ)
  • 既存の方式(折れ線近似)に比較し、パフォーマンス、精度とも大幅に向上

※HPG 2017 の Wolfgang Straßer Award (Best Paper) を受賞した優秀な論文です。今回、高次方程式の解法のみを切り出して参考にしましたが、レイトレーシングやPBRTに関心のある方は全文を一読されることをお奨めします。

準備

Reshetovさんの論文から、ブダン・フーリエの定理、ヴァンサンの定理、リダーズ法やミュラー法がカギとなることはわかりましたが、これらについての日本語の解説が見あたらず、ひとまず勉強がてら、以下を英語版Wikipediaから翻訳して作成しました。

※ ブダンの定理は、内容的には『ブダン・フーリエの定理』とするほうが適切なように思われましたが、今回は英語版からの翻訳であり、原文のとおり『ブダンの定理』としました。
※ Vincentは英語読みだと『ヴィンセント』ですが、フランスの数学者なので他の例にならいフランス語読みの『ヴァンサン』にしました。

実験した方式

対象の9次関数($F_0(u)$ または $F_1(t)$)を $p(x)$ とします。

ステップ1

ブダン・フーリエの定理を用いて、実根の有無を検定する。
$x=0$ におけるフーリエ数列 $\lbrace p(0), p^{(1)}(0), p^{(2)}(0), \dots, p^{(9)}(0) \rbrace$ の符号変化の数と、$x=1$ におけるフーリエ数列 $\lbrace p(1), p^{(1)}(1), p^{(2)}(1), \dots, p^{(9)}(1) \rbrace$ の符号変化の数の差を求める。

  • 0であれば実根なし →処理終了
  • 1であれば実根が1個ある →ステップ3へ
  • 2以上であれば実根が1個以上ある →ステップ2へ

ステップ2

VAG法(Vincent–Alesina–Galuzzi:ヴァンサンの定理の応用手法)を用いて、実根が1個含まれる区間を求める。
区間 $[a, b]$ について、変換された多項式 $ f(x) = (x+1)^9 p \left (\frac{a+bx}{1+x} \right ) $ の係数列の符号変化の数を求める。(初期値は $a=0, b=1$)

  • 0であれば実根なし →区間 $[a, b]$ は棄却
  • 1であれば実根が1個ある →区間 $[a, b]$ をステップ3で処理
  • 2以上であれば実根が1個以上ある →区間 $[a, b]$ を2分割し、反復する

ステップ3

二分法により、実根が1個含まれる区間を1/8に縮小する。
※このステップがないと、次のステップで収束性が悪くなったり、ニュートン・ラフソン法では別の根に収束してしまう等の問題が発生したため、追加しました。1/8は『とりあえず』です。

ステップ4

ニュートン・ラフソン法およびリダーズ法を適用し、それぞれSage Mathで求めた根(真値とする)に収束するか確認する。また、両者の収束性を比較する。

実験方法

VAG法の多項式の係数はSage Mathで求め、その他の計算はExcel上で行い人力でまとめました。(論文にはフーリエ数列を多項式の係数列に変換する方法が提示されていますが、今回は見送りました)

例として、VAG法の多項式の係数を求めるためのSage Mathのコードは以下の通りです。fを求める際、項ごとにいちいち「*(u+1)^9」しているのは、こうしないとsimplify_fullで約分されない場合があるからです。

VAG法の多項式の係数を求める(座標は例①用)
p0x=7; p1x=23; p2x=1; p3x=15
p0y=8; p1y=20; p2y=1; p3y=11
p4x=10; p5x=22; p6x=2; p7x=12
p4y=11; p5y=5; p6y=20; p7y=7
var('t u')
a0=-p0x+3*p1x-3*p2x+p3x
a1=3*p0x-6*p1x+3*p2x
a2=-3*p0x+3*p1x
a3=p0x
b0=-p0y+3*p1y-3*p2y+p3y
b1=3*p0y-6*p1y+3*p2y
b2=-3*p0y+3*p1y
b3=p0y
c0=-p4x+3*p5x-3*p6x+p7x
c1=3*p4x-6*p5x+3*p6x
c2=-3*p4x+3*p5x
c3=p4x
d0=-p4y+3*p5y-3*p6y+p7y
d1=3*p4y-6*p5y+3*p6y
d2=-3*p4y+3*p5y
d3=p4y
x=c0*u^3+c1*u^2+c2*u+c3
y=d0*u^3+d1*u^2+d2*u+d3
p=a3^3*b0^3-a2*a3^2*b0^2*b1+a1*a3^2*b0*b1^2-a0*a3^2*b1^3+a2^2*a3*b0^2*b2-2*a1*a3^2*b0^2*b2-a1*a2*a3*b0*b1*b2+3*a0*a3^2*b0*b1*b2+a0*a2*a3*b1^2*b2+a1^2*a3*b0*b2^2-2*a0*a2*a3*b0*b2^2-a0*a1*a3*b1*b2^2+a0^2*a3*b2^3-a2^3*b0^2*b3+3*a1*a2*a3*b0^2*b3-3*a0*a3^2*b0^2*b3+a1*a2^2*b0*b1*b3-2*a1^2*a3*b0*b1*b3-a0*a2*a3*b0*b1*b3-a0*a2^2*b1^2*b3+2*a0*a1*a3*b1^2*b3-a1^2*a2*b0*b2*b3+2*a0*a2^2*b0*b2*b3+a0*a1*a3*b0*b2*b3+a0*a1*a2*b1*b2*b3-3*a0^2*a3*b1*b2*b3-a0^2*a2*b2^2*b3+a1^3*b0*b3^2-3*a0*a1*a2*b0*b3^2+3*a0^2*a3*b0*b3^2-a0*a1^2*b1*b3^2+2*a0^2*a2*b1*b3^2+a0^2*a1*b2*b3^2-a0^3*b3^3+(-3*a3^2*b0^3+2*a2*a3*b0^2*b1-2*a1*a3*b0*b1^2+2*a0*a3*b1^3-a2^2*b0^2*b2+4*a1*a3*b0^2*b2+a1*a2*b0*b1*b2-6*a0*a3*b0*b1*b2-a0*a2*b1^2*b2-a1^2*b0*b2^2+2*a0*a2*b0*b2^2+a0*a1*b1*b2^2-a0^2*b2^3-3*a1*a2*b0^2*b3+6*a0*a3*b0^2*b3+2*a1^2*b0*b1*b3+a0*a2*b0*b1*b3-2*a0*a1*b1^2*b3-a0*a1*b0*b2*b3+3*a0^2*b1*b2*b3-3*a0^2*b0*b3^2)*x+(a2^3*b0^2-3*a1*a2*a3*b0^2+3*a0*a3^2*b0^2-a1*a2^2*b0*b1+2*a1^2*a3*b0*b1+a0*a2*a3*b0*b1+a0*a2^2*b1^2-2*a0*a1*a3*b1^2+a1^2*a2*b0*b2-2*a0*a2^2*b0*b2-a0*a1*a3*b0*b2-a0*a1*a2*b1*b2+3*a0^2*a3*b1*b2+a0^2*a2*b2^2-2*a1^3*b0*b3+6*a0*a1*a2*b0*b3-6*a0^2*a3*b0*b3+2*a0*a1^2*b1*b3-4*a0^2*a2*b1*b3-2*a0^2*a1*b2*b3+3*a0^3*b3^2)*y+(3*a3*b0^3-a2*b0^2*b1+a1*b0*b1^2-a0*b1^3-2*a1*b0^2*b2+3*a0*b0*b1*b2-3*a0*b0^2*b3)*x^2+(3*a1*a2*b0^2-6*a0*a3*b0^2-2*a1^2*b0*b1-a0*a2*b0*b1+2*a0*a1*b1^2+a0*a1*b0*b2-3*a0^2*b1*b2+6*a0^2*b0*b3)*x*y+(a1^3*b0-3*a0*a1*a2*b0+3*a0^2*a3*b0-a0*a1^2*b1+2*a0^2*a2*b1+a0^2*a1*b2-3*a0^3*b3)*y^2-b0^3*x^3+3*a0*b0^2*x^2*y-3*a0^2*b0*x*y^2+a0^3*y^3
e=p.coefficients(u,sparse=False)
a=0
b=1
u1=(a+b*u)/(1+u)
f=e[9]*u1^9*(u+1)^9+e[8]*u1^8*(u+1)^9+e[7]*u1^7*(u+1)^9+e[6]*u1^6*(u+1)^9+e[5]*u1^5*(u+1)^9+e[4]*u1^4*(u+1)^9+e[3]*u1^3*(u+1)^9+e[2]*u1^2*(u+1)^9+e[1]*u1^1*(u+1)^9+e[0]*(u+1)^9
f=f.simplify_full()
f.coefficients(u,sparse=False)

実験結果

例①(前回の例①と同一)

 inter_ex1.png

曲線0(オレンジ色)の9次方程式を $F_0(u)=0$ とします。
曲線1(青色)の9次方程式を $F_1(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)$の根
1 0.0192470307676812 0.0705390793450691
2 0.0549916382769417 0.0972454495384701
3 0.122057002218322 0.151317532883142
4 0.320989106683592 0.423317135673883
5 0.442530534075141 0.504784914630604
6 0.488155778702660 0.616234378772213
7 0.900475882753311 0.858627731561791
8 0.944674407179144 0.947328062423360
9 0.963072466334768 0.966598091975171

ステップ1

$F_0(u)$ と $F_1(t)$ についてのフーリエ数列は以下の通りです。
両者とも、符号変化の数の差が9であり、実根が1~9の奇数個あることがわかります。
ステップ2に進みます。

$u=0$ $u=1$
$F_0(u)$ 2.91E+06 -1.27E+07
$F_0^{(1)}(u)$ -2.59E+08 -8.09E+08
$F_0^{(2)}(u)$ 1.39E+10 -3.66E+10
$F_0^{(3)}(u)$ -4.81E+11 -1.11E+12
$F_0^{(4)}(u)$ 1.11E+13 -2.20E+13
$F_0^{(5)}(u)$ -1.78E+14 -3.01E+14
$F_0^{(6)}(u)$ 1.99E+15 -2.88E+15
$F_0^{(7)}(u)$ -1.49E+16 -1.88E+16
$F_0^{(8)}(u)$ 6.80E+16 -7.58E+16
$F_0^{(9)}(u)$ -1.44E+17 -1.44E+17
$t=0$ $t=1$
$F_1(t)$ -4.26E+07 7.70E+06
$F_1^{(1)}(t)$ 1.72E+09 5.06E+08
$F_1^{(2)}(t)$ -5.36E+10 2.30E+10
$F_1^{(3)}(t)$ 1.26E+12 6.99E+11
$F_1^{(4)}(t)$ -2.22E+13 1.45E+13
$F_1^{(5)}(t)$ 2.90E+14 2.14E+14
$F_1^{(6)}(t)$ -2.75E+15 2.23E+15
$F_1^{(7)}(t)$ 1.81E+16 1.59E+16
$F_1^{(8)}(t)$ -7.41E+16 6.98E+16
$F_1^{(9)}(t)$ 1.44E+17 1.44E+17

ステップ2

$F_0(u)$ について、VAG法により多項式 $f(x)$ の係数の符号変化の数を求め、区間分割することを反復して実根分離したチャートは以下の通りです。
赤枠は実根がないので棄却、青枠は実根1個で分離完了した区間です。
VAG1.png
$F_1(t)$ については以下の通りです。
VAG2.png

ステップ3

$F_0(u)$ について、分離された区間 $[a,b]$ と、二分法で縮小した区間は以下の通りです。

$a$ $b$ 縮小後$a$ 縮小後$b$
1 0 0.03125 0.015625 0.01953125
2 0.03125 0.0625 0.0546875 0.05859375
3 0.0625 0.125 0.1171875 0.125
4 0.25 0.375 0.3125 0.328125
5 0.4375 0.46875 0.44140625 0.4453125
6 0.46875 0.5 0.484375 0.48828125
7 0.875 0.9375 0.8984375 0.90625
8 0.9375 0.953125 0.943359375 0.9453125
9 0.953125 0.96875 0.962890625 0.96484375

$F_1(t)$ について、分離された区間 $[a,b]$ と、二分法で縮小した区間は以下の通りです。

$a$ $b$ 縮小後$a$ 縮小後$b$
1 0.0625 0.09375 0.0703125 0.07421875
2 0.09375 0.125 0.09375 0.09765625
3 0.125 0.25 0.140625 0.15625
4 0.25 0.5 0.40625 0.4375
5 0.5 0.5625 0.5 0.5078125
6 0.5625 0.625 0.609375 0.6171875
7 0.75 0.875 0.84375 0.859375
8 0.9375 0.953125 0.947265625 0.94921875
9 0.953125 0.96875 0.96484375 0.966796875

ステップ4

$F_0(u)$ についての収束グラフは以下の通りです。(左:ニュートン法、右:リダーズ法)
横軸が反復回数で0は初期値、縦軸は真値(Sage Mathで求めた根)との差の絶対値です。
おしなべてリダーズ法のほうが収束が速いです。ニュートン法はばらつきがあります。
newton_ex1u.png ridders_ex1u.png hanrei.png

$F_1(t)$ についての収束グラフは以下の通りです。
newton_ex1t.png ridders_ex1t.png hanrei.png

例②(前回の例②と同一)

 inter_ex3.png

曲線0(オレンジ色)の9次方程式を $F_0(u)=0$ とします。
曲線1(青色)の9次方程式を $F_1(t)=0$ とします。

$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
$F_0(u)$の根 $F_1(t)$の根
1 -17.4455311309653 -11.7491222197159
2 0.559001039805706 -1.13495840453138
3 1.58591522923000 0.427184956901181
4 4.61782798597658 1.82778197153885
5 5.80036113872243 2.78366617296693
6 -0.584112384659363 - 0.893083467241210*I -10.1061675932554 - 1.15206211156715*I
7 -0.584112384659363 + 0.893083467241210*I -10.1061675932554 + 1.15206211156715*I
8 12.0253252532746 - 13.4619234074658*I 0.528891354675582 - 0.942411682190976*I
9 12.0253252532746 + 13.4619234074658*I 0.528891354675582 + 0.942411682190976*I

ステップ1

$F_0(u)$ と $F_1(t)$ についてのフーリエ数列は以下の通りです。
両者とも、符号変化の数の差が3であり、実根が1または3個あることがわかります。
ステップ2に進みます。

$u=0$ $u=1$
$F_0(u)$ -1.54E+05 8.29E+04
$F_0^{(1)}(u)$ 2.77E+05 8.42E+04
$F_0^{(2)}(u)$ -4.05E+04 -6.18E+05
$F_0^{(3)}(u)$ 7.69E+05 -9.51E+05
$F_0^{(4)}(u)$ -5.03E+06 8.15E+05
$F_0^{(5)}(u)$ 8.16E+06 3.55E+06
$F_0^{(6)}(u)$ -4.67E+06 -4.46E+06
$F_0^{(7)}(u)$ -9.07E+04 4.54E+05
$F_0^{(8)}(u)$ 7.26E+05 3.63E+05
$F_0^{(9)}(u)$ -3.63E+05 -3.63E+05
$t=0$ $t=1$
$F_1(t)$ -3.50E+03 3.19E+03
$F_1^{(1)}(t)$ 1.05E+04 4.94E+03
$F_1^{(2)}(t)$ -1.10E+04 -5.02E+03
$F_1^{(3)}(t)$ -1.37E+04 -3.56E+04
$F_1^{(4)}(t)$ 1.87E+05 -1.54E+05
$F_1^{(5)}(t)$ -4.54E+05 2.92E+04
$F_1^{(6)}(t)$ -1.90E+05 1.37E+06
$F_1^{(7)}(t)$ 9.53E+05 2.22E+06
$F_1^{(8)}(t)$ 1.09E+06 1.45E+06
$F_1^{(9)}(t)$ 3.63E+05 3.63E+05

ステップ2

$F_0(u)$ と $F_1(t)$ について、VAG法により多項式 $f(x)$ の係数と、符号変化の数を求めた結果は以下の通りです。
両者とも、符号変化の数が1であり、区間分割・反復せずに実根分離完了となります。

$F_0(u)$の$f(x)$の係数 $F_1(t)$の$f(x)$の係数
0 8.29E+04 3.19E+03
1 6.62E+05 2.37E+04
2 2.00E+06 7.27E+04
3 2.60E+06 1.18E+05
4 2.25E+05 1.01E+05
5 -3.75E+06 2.44E+04
6 -5.18E+06 -4.15E+04
7 -3.34E+06 -4.77E+04
8 -1.11E+06 -2.10E+04
9 -1.54E+05 -3.50E+03

ステップ3

$F_0(u)$ について、分離された区間 $[a,b]$ と、二分法で縮小した区間は以下の通りです。

$a$ $b$ 縮小後$a$ 縮小後$b$
0 1 0.5 0.625

$F_1(t)$ について、分離された区間 $[a,b]$ と、二分法で縮小した区間は以下の通りです。

$a$ $b$ 縮小後$a$ 縮小後$b$
0 1 0.375 0.5

ステップ4

$F_0(u)$ についての収束グラフは以下の通りです。(左:ニュートン法、右:リダーズ法)
newton_ex3u.png ridders_ex3u.png

$F_1(t)$ についての収束グラフは以下の通りです。
newton_ex3t.png ridders_ex3t.png

例③(前回DKA法で敗北した例③と同一)

 inter_ex5.png

曲線0(オレンジ色)の9次方程式を $F_0(u)=0$ とします。
曲線1(青色)の9次方程式を $F_1(t)=0$ とします。

$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)$の根
1 -567115.622411213 -248197.018538645
2 0.557665338289116 -1.03330670884847
3 1.61807431932095 0.446194297282387
4 5.16799108368957 1.86453498958389
5 5.79814264427882 2.90335369440187
6 -0.570941565478316 - 0.800872586604778*I -213844.115726888 - 446351.819141408*I
7 -0.570941565478316 + 0.800872586604778*I -213844.115726888 + 446351.819141408*I
8 283559.811210479 - 491136.142417216*I 0.534608073955102 - 0.882297985704637*I
9 283559.811210479 + 491136.142417216*I 0.534608073955102 + 0.882297985704637*I

ステップ1

$F_0(u)$ と $F_1(t)$ についてのフーリエ数列は以下の通りです。
$F_0(u)$ については符号変化の数の差が1であり、実根が1個あることがわかるのでステップ3に進みます。
$F_1(t)$ については符号変化の数の差が3であり、実根が1または3個あることがわるのでステップ2に進みます。

$u=0$ $u=1$
$F_0(u)$ -1.36E+05 8.81E+04
$F_0^{(1)}(u)$ 2.17E+05 1.06E+05
$F_0^{(2)}(u)$ 5.98E+04 -5.95E+05
$F_0^{(3)}(u)$ 7.87E+05 -1.16E+06
$F_0^{(4)}(u)$ -5.06E+06 5.37E+05
$F_0^{(5)}(u)$ 7.46E+06 3.73E+06
$F_0^{(6)}(u)$ -3.73E+06 -3.73E+06
$F_0^{(7)}(u)$ 2.18E-04 2.18E-04
$F_0^{(8)}(u)$ 1.83E-08 8.02E-09
$F_0^{(9)}(u)$ -1.03E-08 -1.03E-08
$t=0$ $t=1$
$F_1(t)$ -4.59E+03 3.19E+03
$F_1^{(1)}(t)$ 1.45E+04 4.94E+03
$F_1^{(2)}(t)$ -2.06E+04 -5.02E+03
$F_1^{(3)}(t)$ -1.43E+04 -3.24E+04
$F_1^{(4)}(t)$ 3.19E+05 -1.48E+05
$F_1^{(5)}(t)$ -1.09E+06 1.56E+05
$F_1^{(6)}(t)$ 1.24E+06 1.24E+06
$F_1^{(7)}(t)$ 5.03E+01 5.03E+01
$F_1^{(8)}(t)$ 7.75E-04 7.75E-04
$F_1^{(9)}(t)$ 7.74E-09 7.74E-09

ステップ2

$F_1(t)$ について、VAG法により多項式 $f(x)$ の係数と、符号変化の数を求めた結果は以下の通りです。
符号変化の数が1であり、区間分割・反復せずに実根分離完了となります。

$F_1(t)$の$f(x)$の係数
0 3.19E+03
1 2.37E+04
2 7.27E+04
3 1.17E+05
4 9.83E+04
5 1.66E+04
6 -5.40E+04
7 -5.95E+04
8 -2.68E+04
9 -4.59E+03

ステップ3

$F_0(u)$ について、分離された区間 $[a,b]$ と、二分法で縮小した区間は以下の通りです。

$a$ $b$ 縮小後$a$ 縮小後$b$
0 1 0.5 0.625

$F_1(t)$ について、分離された区間 $[a,b]$ と、二分法で縮小した区間は以下の通りです。

$a$ $b$ 縮小後$a$ 縮小後$b$
0 1 0.375 0.5

ステップ4

$F_0(u)$ についての収束グラフは以下の通りです。(左:ニュートン法、右:リダーズ法)
newton_ex6u.png ridders_ex6u.png

$F_1(t)$ についての収束グラフは以下の通りです。
newton_ex6t.png ridders_ex6t.png

例④(交点をカーブセグメント端部に寄せた)

 inter_ex6.png

曲線0(オレンジ色)の9次方程式を $F_0(u)=0$ とします。
曲線1(青色)の9次方程式を $F_1(t)=0$ とします。

$F_0(u)$の係数 $F_1(t)$の係数
0 -125000 125000
1 1237500 -967500
2 -5681250 3396150
3 15354354 -6779814
4 -26309835 8094645
5 28696293 -5454639
6 -19066065 1614507
7 6774552 45891
8 -988191 -2232
9 106651 -107
$F_0(u)$の根 $F_1(t)$の根
1 0.996305429060318 0.0437868412108868
2 0.0752756536928661 - 0.147597702909089*I -0.0327786776634256 - 0.0179621297154259*I
3 0.0752756536928661 + 0.147597702909089*I -0.0327786776634256 + 0.0179621297154259*I
4 1.00075116412009 - 0.775815852180924*I 1.14224445195882 - 0.737090574225273*I
5 1.00075116412009 + 0.775815852180924*I 1.14224445195882 + 0.737090574225273*I
6 1.52882423470616 - 1.42568576467889*I 1.29071708183460 - 1.04366047514815*I
7 1.52882423470616 + 1.42568576467889*I 1.29071708183460 + 1.04366047514815*I
8 1.84699623295072 - 1.02025730953808*I 1.44792372326456 - 0.807302190695345*I
9 1.84699623295072 + 1.02025730953808*I 1.44792372326456 + 0.807302190695345*I

ステップ1

$F_0(u)$ と $F_1(t)$ についてのフーリエ数列は以下の通りです。
$F_0(u)$ については符号変化の数の差が5であり、実根が1~5の奇数個あることがわかります。
$F_1(t)$ については符号変化の数の差が7であり、実根が1~7の奇数個あることがわかります。
ステップ2に進みます。

$u=0$ $u=1$
$F_0(u)$ 1.07E+05 -9.91E+02
$F_0^{(1)}(u)$ -9.88E+05 -2.69E+05
$F_0^{(2)}(u)$ 1.35E+07 -3.70E+05
$F_0^{(3)}(u)$ -1.14E+08 -2.02E+06
$F_0^{(4)}(u)$ 6.89E+08 -1.22E+07
$F_0^{(5)}(u)$ -3.16E+09 7.20E+06
$F_0^{(6)}(u)$ 1.11E+10 -1.90E+08
$F_0^{(7)}(u)$ -2.86E+10 -1.42E+09
$F_0^{(8)}(u)$ 4.99E+10 4.54E+09
$F_0^{(9)}(u)$ -4.54E+10 -4.54E+10
$t=0$ $t=1$
$F_1(t)$ -1.07E+02 7.19E+04
$F_1^{(1)}(t)$ -2.23E+03 6.69E+04
$F_1^{(2)}(t)$ 9.18E+04 2.80E+05
$F_1^{(3)}(t)$ 9.69E+06 1.99E+06
$F_1^{(4)}(t)$ -1.31E+08 5.08E+06
$F_1^{(5)}(t)$ 9.71E+08 3.66E+07
$F_1^{(6)}(t)$ -4.88E+09 2.90E+08
$F_1^{(7)}(t)$ 1.71E+10 7.87E+08
$F_1^{(8)}(t)$ -3.90E+10 6.35E+09
$F_1^{(9)}(t)$ 4.54E+10 4.54E+10

ステップ2

$F_0(u)$ について、VAG法により多項式 $f(x)$ の係数の符号変化の数を求め、区間分割することを反復して実根分離したチャートは以下の通りです。
赤枠は実根がないので棄却、青枠は実根1個で分離完了した区間です。
VAG3.png

$F_1(t)$ について、VAG法により多項式 $f(x)$ の係数と、符号変化の数を求めた結果は以下の通りです。符号変化の数が1であり、区間分割・反復せずに実根分離完了となります。

$F_1(t)$の$f(x)$の係数
0 7.19E+04
1 5.80E+05
2 2.19E+06
3 4.81E+06
4 6.48E+06
5 5.06E+06
6 1.86E+06
7 2.42E+04
8 -3.20E+03
9 -1.07E+02

ステップ3

$F_0(u)$ について、分離された区間 $[a,b]$ と、二分法で縮小した区間は以下の通りです。

$a$ $b$ 縮小後$a$ 縮小後$b$
0.5 1 0.9375 1

$F_1(t)$ について、分離された区間 $[a,b]$ と、二分法で縮小した区間は以下の通りです。

$a$ $b$ 縮小後$a$ 縮小後$b$
0 1 0 0.125

ステップ4

$F_0(u)$ についての収束グラフは以下の通りです。(左:ニュートン法、右:リダーズ法)
newton_ex5u.png ridders_ex5u.png

$F_1(t)$ についての収束グラフは以下の通りです。
newton_ex5t.png ridders_ex5t.png

$F_1(t)$ についての収束性が他より悪いのは、対象の関数・区間内に極小点が含まれることが影響したようです。(以下のグラフ参照。赤丸:根、青三角:極小点)
対象範囲を絞る方法の改善が必要と思われます。
graph_ex5t.png

課題

  • 対象範囲を絞る方法を検討する(ベジェクリッピングを1巡だけ適用、なども良いかも)
  • 他の求根法(ミュラー法など)も試してみる
  • VAG法の多項式の係数列を、フーリエ数列から変換する方法を試してみる
  • 複数の根が近接、あるいは重なっている(重根)場合の対応をきちんと設計する
0
0
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
0
0

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?