FCMP.SOLVEで方程式の実数解をガウス・ニュートン法で求める

 まずSASの公式に乗ってる例から

y = 1 / √x

でy=20となるxの値を求めてみます

グラフにすると

なので，まあ　だいたい値は0.0023～26あたりかなぁって感じですが

proc fcmp outlib=work.common.solve ;

/*解きたい式を定義*/

   function f1(x);

     y = 1/sqrt(x);

      return(y);

   endsub;

 

/*yを指定して、xを求めて戻すルーチン*/

  subroutine solver(y,x);

   outargs x  ;

   x=solve("f1",{.}, y , .);

  endsub ;

run;

options cmplib=work.common;

data wk1;

 call missing(x);

 call solver(20,x);

run;

proc print data=wk1;

run;

👆
と，このように，あらかじめ解きたい式をFCMPで指定し，
FCMP.SOLVE関数の第一引数に　解きたい関数名，第三引数に目標値（ここでいうところの左辺，y）を指定します，
第二引数には計算のための5つのパラメータを配列で与えるのですが，{.}とするとデフォルト設定になります
最後の引数は，第一引数にとった関数に追加引数をここで与えることができます

つぎに3次方程式

を解いてみます

グラフにすると

👆どうも実数解は3つありそうです．ぱっとみ，交点は1,2,3らへん

この場合，初期値の設定によって，収束する解が異なるので，コールルーチンをちょっと改造して，パラメータを細かく指定できるようにして実行してみます

proc fcmp outlib=work.common.solve ;

/*解きたい式を定義*/

   function f1(x);

      return(x**3 -6*x**2 + 11*x -6);

   endsub;

 

/*yを指定して、xを求めて戻すルーチン*/

  subroutine solver(y,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

  /*設定パラメータ */

   array opts[5] initial_value absolute_criterion relative_criterion maximum_iterations solve_status;

   outargs x ,solve_status,y ;

   x=solve("f1"

               ,opts

               , 0

               , .);

   put x= solve_status=;

  endsub ;

run;

options cmplib=work.common;

data wk1;

 call missing(of x);

 absolute_criterion=1.0e-12;

 relative_criterion=1.0e-6;

 maximum_iterations=100;

 solve_status=.;

 y=0;

 initial_value=0.001;

 call solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

 output;

 initial_value=1.9;

 call solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

 output;

 initial_value=2.9;

 call solver(20,x,initial_value,absolute_criterion,relative_criterion,maximum_iterations,solve_status);

 output;

run;

proc print data=wk1;

run;

グラフから，おそらくｘ＝1,2,3あたりに解がありそうなので
初期値(initial_value)を，0.001，1.9, 2.9として回してみると

👆
やはり ｘ=1,2,3が解でした．solve_status=0は収束していることを示します