頭と尻尾はくれてやる!

パソコンおやじのiPhoneアプリ・サイト作成・運営日記


ニュートン法で三次方程式を解く

最近ベジェ曲線と戯れてるんだ。Bezier曲線ね。でもすぐに忘れてベジェかペジェかどっちだっけ?ってなるんだけどさ。
iOSでベジェ曲線をいじるのは初めてなので試行錯誤しながら遊んでいるんだ。
やりたいのが3次のベジェ曲線を描いて、途中の点の座標を得たいんだよ。

そのあたりで三次方程式を解く必要が出てきてさ。二次方程式だと公式あるじゃない?中学か高校で覚えたでしょ?あんなのが三次のでもあったんじゃなかったかなって思って調べるとさ、こんなページがあったんだ。

三次方程式の解の公式 [物理のかぎしっぽ]

あー、あるじゃない。カルダノの公式って言うんだね、聞いたことないや。ともかくさっそくやってみたんだけど、いくらやっても答えが得られないんだ。
よくよく調べると計算の途中でルートの中が負になることがあるんだ。虚数が出てくるんだよ。虚数なんて扱えないからそれ以降はNaNになっちゃうんだ。
面白いことに、答えが実数でも計算の途中で虚数が出てくるとかなんとか。不思議だね。これをコードで実装する(実部を虚部をわけてうんぬん、、、とかになるのかな?)のはちょっと難儀するねえってことで結局繰り返し計算するニュートン法ってのをやってみたよ。

x^3 + ax^2 + bx + c = 0
↑こういう三次方程式があるとするよ。これを解きたいってことは
y = x^3 + ax^2 + bx + c
って式のyがゼロになる時のxを求めるってことだよね。

ニュートン法概略

↑ざっくりとした図だけど、こんなのだ。
初期値x0を与えてその時のyを求める。これがゼロならいいわけだけど、ゼロじゃない(というか設定したしきい値より大きい場合)ならそのxでの接線がx軸と交わるところを次のxとして、収束するまで繰り返し計算するわけだ。

コードはこんな感じ。xじゃなくてtになってるけど。
{
    double limitDt = 0.001;//しきい値
    double t = t0;//初期値
    
    while (1) {
        double dt = pow(t, 3) + a*pow(t, 2) + b*t+c;
        
        if (ABS(dt) < limitDt) break;//判定
        
        //tにおける傾き
        double slope = 3.f*pow(t,2 ) + 2.f*a*t+b;
        
        //新しいtを得る
        t = (slope*t - dt )/slope;
        
    }
}
実際には繰り返し回数が無限になったりしないようにチェックとかしてるけど、基本は上のような感じで。
初期値の与え方やしきい値の設定にもよるけど、俺コードの場合だとだいたい3回程度で収束してる。
カルダノの公式を頑張って解くより多分楽チンだと思うなあ。

ベジェ曲線

↑ってことで、無事にベジェ曲線の座標を得ることができそうだよ。
これはテスト中の画面だけど、三次のベジェ曲線を描いて、タップしたところのx座標を拾ってグレーの縦線を描いて、その縦線とベジェ曲線の交点にピンクの点を描いてるんだけど、x座標からt(0≦t≦1)を求める時に三次方程式を解いてるのよ。

<< Google Driveのアップデート版がアルバムになる!  TopPage  アプリが落ちるの"落ちる"は英語でなんて言う? >>

コメント


管理者にだけ表示を許可する
 

トラックバック

トラックバックURL
http://ringsbell.blog117.fc2.com/tb.php/803-f779f811




Copyright ©頭と尻尾はくれてやる!. Powered by FC2 Blog. Template by eriraha.

FC2Ad