初心者によるMATLABメモ

検索エンジンから来た方は読みたい記事が表示されていない可能性があります。
その場合、左メニューのフォームでブログ内を検索すると見たい記事を見れると思います。
もしくは、カテゴリー名をクリックしていただくと各カテゴリーの記事のリストが見れます。

上記の広告は１ヶ月以上更新のないブログに表示されています。
新しい記事を書く事で広告が消せます。

たまたまここでで見たんですが、↓だそうです。

2月最後の授業で言ったように､Matlabのoptimization toolboxを持っている場合､この中の関数fsolveを用いることによって非線形の連立方程式の解を求めることができます。しかし、その後調べたところ､数値計算に関する教科書を著しているKenn Juddは彼のホームページでこの関数はイマイチ（”poor”）であると述べ(リンク先参照)、それよりもChristopher Simsの書いた関数コードであるcsolve.mの利用を推奨しています。…
…
参考：コメント：Matlab関数fsolveについて

で、実際にリンクをたどってみると、確かに書いてありますね。そんなものを売るな！という話もありますが、まぁそれはそれ。無料の方が強力なら、それもよしです。

A Matlab Nonlinear Equation Solver The nonlinear equation solver in Matlab, fsolve, is a poor one based on minimizing the sum of squares of the functions. Chris Sims offers an improvement; see his ftp page.…
…
参考：Chapter 5 public code

とりあえずcsolve()が欲しいなら、ここからダウンロードできます。

以下で、csolve()の使い方を書いておきます。
非線形の方程式ソルバーcolve()は，準ニュートン法ベースで、通常の解法で極値，振動などで解が求まらない場合，ランダム方向に解を検索することで解を改善しようとします．
実際の使い方は下記です。

[x,rc] = csolve(FUN,x0,gradfun,crit,itmax)
[x,rc] = csolve(FUN,x0,gradfun,crit,itmax,varargin)

% [引数]
FUN        :    解かれる非線形システム方程式。
            FUNはベクトルxを引数として計算した非線形方程式のベクトルFを返します。
            ベクトルFはベクトルxと常に同じサイズである必要があります。
            関数FUNは、関数ハンドルを使って指定することができます。
x0        :    解の探索のための初期値．
gradfun    :    ヤコビ行列の関数．
crit    :    FUNが返す絶対値の合計がこの値より小さいならば、解が収束したと判断されます。
itmax    :    許可する繰り返しの最大回数．もし繰り返しの回数がこの値に達した場合、第2出力rcに4が返ります。
varargin:    関数FUNとgradfunに渡す、ベクトルx以外のパラメータを指定します。

この関数は実際に使った例を示します。
ここでは、下記の非線形方程式を解いてみます。

この式はfsolveのときと同じ方程式です。
方程式の関数は下記です。

% 方程式の関数
function F = funcs(x) 
    F = zeros(numel(x),1);
    F(1) = x(1)^2+x(2)^2-4;
    F(2) = x(1)*x(2)-x(3);
    F(3) = x(1)*x(3)-1;
    F = F(:);

% 関数funcsのヤコビ行列
function J = funcsj(x) 
    J = zeros(numel(x));
    J(1,1) = 2*x(1);
    J(2,1) = x(2);
    J(3,1) = x(3);

    J(1,2) = 2*x(2);
    J(2,2) = x(1);
    J(3,2) = 0;

    J(1,3) = 0;
    J(2,3) = -1;
    J(3,3) = x(1);

実行例は、下記です。fsolveのときと同じ解になっていることがわかります。

>> [x,exitflag] = mcsolve(@funcs, ones(3,1), @funcsj, 1e-12, 100, 4, 1)
itct 1, af 1.0208, lambda 0.36, rc 0
   x        0.64       1.72       1.36 
   f      -0.632    -0.2592    -0.1296 
itct 2, af 0.0452359, lambda 1, rc 0
   x    0.735407    1.86822    1.35976 
   f   0.0310718  0.0141413 -2.28221e-005 
itct 3, af 8.75519e-005, lambda 1, rc 0
   x    0.733074    1.86082     1.3641 
   f  6.0167e-005 1.72542e-005 -1.01306e-005 
itct 4, af 3.64242e-010, lambda 1, rc 0
   x    0.733077    1.86081    1.36411 
   f  2.9954e-010 -4.20057e-011 2.26965e-011 
itct 5, af 8.88178e-016, lambda 1, rc 0
   x    0.733077    1.86081    1.36411 
   f  8.88178e-016          0          0 
x =
    0.7331
    1.8608
    1.3641


exitflag =
     0

ただし、この関数csolve()は少々使いづらい部分もありますね。そのあたりは次のエントリに書きたいと思います。

参考：コメント：Matlab関数fsolveについて
参考：csolve()のダウンロード先
参考：【文献調査】BFGSの基礎
参考：初心者によるMATLABメモ｜ 非線形方程式を解く(fsolve)