sub_set01へ戻る

名称 ・・・・・・・

3次式の解法(カルダノ方式)

ルーチン名 ・・・・・・・

Cardano

スクリプトNo ・・・・・・・ sub_set01- 015
file名 ・・・・・・・ ScriptSet01.lzh
sub_set01のページからセットごとダウンロードしてください。
概略 ・・・・・・・

Cardanoの公式を使い3次式について各項の係数から解を求める。

使用前提 ・・・・・・・

最高次数(3次)の項の係数を1としとした上で呼び出してください。
結果(戻り値)を受ける通常のバリアント変数を用意(宣言)しておく。

依存関係 ・・・・・・・

特に無し

引数説明 ・・・・・・・

a
b
c

(in)2次式の2次の項の係数。
(in)2次式の1次の項の係数。
(in)2次式の定数項。

戻り値説明 ・・・・・・・

3次式の解(0〜3個)戻る。
要素0は判定用に使われるので、要素1以降に解の個数に応じて結果(解)が格納される。

変数説明 ・・・・・・・

D
nd
t()
q,p,u,v,ux

判別式の結果用の変数
解の数を保管する変数
処理の成否(要素0)と解(要素1以降)が格納される。
計算過程で使用される一時変数

注意点 ・・・・・・・

係数の扱いが1次式の解法,2次式の解法と違い、3次の項の係数を1であるののとして処理するので予め各係数を3次項の係数で割ってから呼び出してください。他のものと別に作成したので整合性が無くてすいません、実際に使用するときは引数に3次項の係数を引数として追加し、サブルーチン内で処理するように改良した方が再利用しやすいかと思います。
複数解の可能性があるので変数tの配列要素数は動的に確保する。
解が無い場合も想定して、戻り値の配列要素0には成否結果を格納するようにした。

 
スクリプト
'--------------3次式の解法(カルダノ方式)
Function Cardano(a,b,c)
'3次式の解をカルダノの公式で求める。
'戻値---解(tパラメータ) 配列(1次元)0番目は成否に使用する。
'引数---
' a   :(in)3次式の2次の項の係数
' b   :(in)3次式の1次の項の係数
' c   :(in)3次式の定数項

  Dim D                  '判別式の結果用の変数
  Dim q,p,u,v,ux         '一時変数
  Dim nd                 '解の数を保管する変数
  Dim tt()               '解(tパラメータ)用変数

  p=-((a/3)^2)+(b/3)
  q=((a/3)^3)-(a/3)*(b/2)+(c/2)
  D=round(((q^2)+(p^3))*1000000)/1000000

  if D>0 then        '解が1つの場合
    nd=1
    Redim tt(nd)
    tt(0)=true
    if -q+sqr(q^2+p^3)>=0 then
      u=(-q+sqr(q^2+p^3))^(1/3)
    else
      u=-((-(-q+sqr(q^2+p^3)))^(1/3))
    end if
    if -q-sqr(q^2+p^3)>=0 then
      v=(-q-sqr(q^2+p^3))^(1/3)
    else
      v=-((-(-q-sqr(q^2+p^3)))^(1/3))
    end if
    tt(1)=u+v-a/3
  elseif D=0 then    '解が2つの場合
    nd=2
    Redim tt(nd)
    tt(0)=true
    if (-q)>=0 then
      tt(1)=2*(-q)^(1/3)-a/3
      tt(2)=-((-q)^(1/3))-a/3
    else
      tt(1)=-2*(-q)^(1/3)-a/3
      tt(2)=-(-q)^(1/3)-a/3
    end if
  elseif D<0 then    '解が3つの場合
    nd=3
    Redim tt(nd)
    tt(0)=true
    ux=q/(p*(-p)^(1/2))
    u=XShade.Acos(ux)
    tt(1)=2*sqr(-p)*cos(u/3)-a/3
    tt(2)=2*sqr(-p)*cos((u/3)+(2*3.14159265358979)/3)-a/3
    tt(3)=2*sqr(-p)*cos((u/3)+(4*3.14159265358979)/3)-a/3
  end if

  Cardano=tt
End Function
使用例
主サブルーチンの例題を参照してください。

sub_set01へ戻る