FC2ブログ

初心者によるMATLABメモ

3次元形状やってんのになぜかMATLABを使わされることになった人のメモ。

 

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

[数値解析] 非線形方程式を解く(準ニュートン法:csolve)

たまたまここでで見たんですが、↓だそうです。
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)
スポンサーサイト



テーマ:MATLAB - ジャンル:コンピュータ

このエントリに付けられたタグ|MATLABOptimizationToolbox方程式非線形方程式fsolvecsolve

[Optimization Toolbox] 非線形方程式を解く2 (fsolve)

このエントリは、こちらのエントリからの続きで,関数がx以外の引数を持つ場合や,オプションを指定してOptimization Toolboxのfsolveを用いて,非線形方程式を解く方法です。
先ほども書きましたが、この方法は、Optimization Toolboxを購入している方しかできません。
(持ってない人はfsolveほど強力ではないですがこちらのニュートン・ラプソン法を参照してください。)

使い方は、こちらのエントリを参照してもらうとして,今回は関数がx以外の引数を持つ場合や,オプションを指定して自作のヤコビ行列を用いて計算する方法です。

ここでも、引き続き下記の非線形な連立方程式を解いてみます。


まずは前回と異なり,関数funcs()がx以外にも引数を持つ場合の例です。
これは簡単に書くと,下記のように4つ目移以降の引数として与えてやればいいだけです。
x = fsolve(@fun,x0,options,fun()に渡すほかの引数...)

まず,意味はないですが,関数funcsを下記のように定義します。a=4で上記の例題の方程式になります。
function F = funcs(x,a) 
    if nargin~=2 a=4; end;
    F = zeros(numel(x),1);
    F(1) = x(1)^2+x(2)^2-a;
    F(2) = x(1)*x(2)-x(3);
    F(3) = x(1)*x(3)-1;
    F = F(:);

実際にfsolveの4つ目の引数として,4を与えて,実行してみると、このように前回出た解x=[0.7331, 1.8608, 1.3641]と一致し,ちゃんと引数が関数に渡っていつることが分かると思います。
>> x = fsolve(@funcs, ones(3,1), opt, 4)

Optimization terminated: first-order optimality is less than options.TolFun.

x =
    0.7331
    1.8608
    1.3641

続いて,ここではヤコビアン行列も関数として実装し、そのヤコビ行列を用いて方程式を解いてみます。
ヤコビアン行列の説明は省略させてもらいますが,上の方程式のヤコビ行列Jは下記になります。


fsolveでは,ヤコビアン行列を指定する場合,ここに書いてあるとおり方程式の関数を下記の形式で定義する必要があります。
function [F,J] = myfun(x)
F = ...          % 方程式の関数
if nargout > 1   % 出力引数が2つ要求するときにだけヤコビ行列を計算する
   J = ...   % ヤコビ行列の計算
end

この形式に倣って,関数funcsj()を実装すると下記になります。
function [F,J] = funcsj(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(:);

    if nargout>1 
        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);
    end;

このヤコビアン行列を使用して方程式を解くには,関数をこの形で実装するだけでは不十分で,さらにfsolve()にオプションを指定してやる必要があり、具体的にはoptimset()で最適化オプションパラメータ構造体を弄ってやる必要があります。
以下に結果とまとめて書いてしまいますが,太文字になっているところで,先ほど定義したヤコビアン行列を使用するオプションをセットしています。
これをセットしない場合は,上記の関数funcsj()内のif nargout>1~endのところは計算されません。
他のオプションも気になる人は,こちらを参照してください。
opt = optimset('Display','iter');  % 最適化オプションパラメータ構造体の生成&繰り返し計算ごとに途中経過を表示するように設定
opt = optimset(opt, 'TolX', 1e-12); % 解xに関する許容誤差を小さくし、真の解に近づける。
opt = optimset(opt, 'TolFun', 1e-10);    % 関数値に関する許容誤差を小さくし、真の解に近づける。
opt = optimset(opt, 'Jacobian','on');    % 自分で定義したヤコビアン行列を使用する。
[x,fval,exitflag,output] = fsolve(@funcsj, ones(3,1), opt)
--以下実行結果
                                         Norm of      First-order   Trust-region
 Iteration  Func-count     f(x)          step         optimality    radius
     0          1               4                             4               1
     1          2        0.545056              1           2.58               1
     2          3      0.00135599       0.187596          0.139             2.5
     3          4     7.4367e-009      0.0102639       0.000322             2.5
     4          5    4.41969e-019   2.68786e-005      2.34e-009             2.5
     5          6    7.88861e-031    2.5354e-010      3.31e-015             2.5
Optimization terminated: first-order optimality is less than options.TolFun.

x =
    0.7331
    1.8608
    1.3641

fval =
  1.0e-015 *
    0.8882
         0
         0

exitflag =
     1

output = 
       iterations: 5
        funcCount: 6
        algorithm: 'trust-region dogleg'
    firstorderopt: 3.3055e-015
          message: 'Optimization terminated: first-order optimality is less than options.TolFun.'

また,ここまでの話を見るとさも方程式の関数を別ファイルで書かないといけないようにも見えますが,別にそういったことはなく無名関数と入れ子関数を使っても問題ありません。
むしろ変数を共有させるために,僕は入れ子関数を使っていることが多いです。
ここらへんはまぁ人それぞれ趣味になってくるのでどうでもいいことですが、一応書いておきます。

参考:Optimization Toolbox User's Guide - fsolve
参考:MATLAB Function Reference - optimset
参考:Optimization Toolbox User's Guide - 最適化オプション
参考:Optimization Toolbox User's Guide - 無名関数と入れ子関数によるグローバル変数の回避
参考:初心者によるMATLABメモ| 非線形方程式を解く(fsolve)
このエントリに付けられたタグ|MATLABOptimizationToolbox方程式非線形方程式fsolveoptimsetヤコビ行列ヤコビアン行列

[Optimization Toolbox] 非線形方程式を解く(fsolve)

Optimization Toolboxのfsolveを用いて,非線形方程式を解く方法です。
ちなみに当然ながら、この方法は、Optimization Toolboxを購入している方のみしかできません。
(持ってない人はfsolveほど強力ではないですがこちらのニュートン・ラプソン法を参照してください。)

まずは使い方。Optimization Toolboxのだいたいの関数がこの使い方なので,Optimization Toolboxを使っているうちに自然に使えるようになると思います。
x = fsolve(@fun,x0) … 初期値x0を用いて関数fun()で定義された非線形方程式を解き、解xを得る。
x = fsolve(@fun,x0,options) … 初期値x0を用い、オプションを構造体optionsで指定して関数fun()で定義された非線形方程式を解を解き、解xを得る。
x = fsolve(@fun,x0,options,fun()に渡すほかの引数...)
[x,fval,exitflag,output,jacobian] = fsolve(...)

このfsolveは解に収束したかどうかはメッセージが出るし、よく出来ているので説明は省略させてもらって実際に例題を解いてみます。

ここでは、下記の非線形方程式を解いてみます。


右辺がすべてゼロになるように整理すると,左辺を下記のような列行列であるとみなせます。
これを関数funcs()として定義します。


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(:);

これを適当に初期値をすべて1として与えて,実行して見た例が下記です。
xが解。fvalは解xを与えたときのfuncs()の値でほぼ0に近くなっています。
exitflagについてはここに書いてある。表を参考にしてください。この場合はexitflag=1になっており,関数が解 xに収束したことがわかります。
続いてoutputは,解を求める際の反復処理に関する情報を持つ構造体です。これについてもexitflag同様ここを参照してください。
最後にjacobianには,解xのときのfuncs()のヤコビ行列です。
ここでは,5つの出力をすべて取得していますが,別にx=fsolve()として解xだけ取得しても構いません。
>> [x,fval,exitflag,output,jacobian] = fsolve(@funcs, ones(3,1)) 
Optimization terminated: first-order optimality is less than options.TolFun.
x =
    0.7331
    1.8608
    1.3641

fval =
  1.0e-009 *
    0.6531
   -0.1214
   -0.0406

exitflag =
     1

output = 
       iterations: 4
        funcCount: 20
        algorithm: 'trust-region dogleg'
    firstorderopt: 2.3414e-009
          message: 'Optimization terminated: first-order optimality is less than options.TolFun.'

jacobian =
    1.4662    3.7216         0
    1.8608    0.7331   -1.0000
    1.3641         0    0.7331

ちなみに、fsolve()のデフォルトのオプションは下記になっています。
これについてはすいませんが,あまりに長くなるのでここでは説明はしません。知りたい人は,こちらを参照してください。
>> optimset('fsolve')
ans = 
                   Display: 'final'
               MaxFunEvals: '100*numberofvariables'
                   MaxIter: 400
                    TolFun: 1.0000e-006
                      TolX: 1.0000e-006
               FunValCheck: 'off'
                 OutputFcn: []
                  PlotFcns: []
           ActiveConstrTol: []
            BranchStrategy: []
           DerivativeCheck: 'off'
               Diagnostics: 'off'
             DiffMaxChange: 0.1000
             DiffMinChange: 1.0000e-008
         GoalsExactAchieve: []
                GradConstr: []
                   GradObj: []
                   Hessian: []
                  HessMult: []
               HessPattern: []
                HessUpdate: []
           InitialHessType: []
         InitialHessMatrix: []
                  Jacobian: 'off'
                 JacobMult: []
              JacobPattern: 'sparse(ones(jrows,jcols))'
                LargeScale: 'off'
        LevenbergMarquardt: []
            LineSearchType: 'quadcubic'
                  MaxNodes: []
                MaxPCGIter: 'max(1,floor(numberofvariables/2))'
                MaxRLPIter: []
                MaxSQPIter: []
                   MaxTime: []
             MeritFunction: []
                 MinAbsMax: []
       NodeDisplayInterval: []
        NodeSearchStrategy: []
          NonlEqnAlgorithm: 'dogleg'
        NoStopIfFlatInfeas: []
      PhaseOneTotalScaling: []
            Preconditioner: []
          PrecondBandWidth: 0
            RelLineSrchBnd: []
    RelLineSrchBndDuration: []
          ShowStatusWindow: []
                   Simplex: []
                    TolCon: []
                    TolPCG: 0.1000
                 TolRLPFun: []
               TolXInteger: []
                  TypicalX: 'ones(numberofvariables,1)'

ちなみに,ちょっとまだまだ長くなりそうなので,オプションを変更して非線形方程式を解く方法は,次のエントリに続きたいと思います。
あと、オプションを弄るoptimsetとかの説明も機会があれば書きたいところです。

参考:Optimization Toolbox User's Guide - fsolve
参考:Optimization Toolbox User's Guide - 最適化オプション
参考:初心者によるMATLABメモ| 非線形方程式を解く2 (fsolve)
参考:初心者によるMATLABメモ| 非線形方程式を解く(ニュートン・ラプソン法)
このエントリに付けられたタグ|MATLABOptimizationToolbox方程式非線形方程式fsolve

[Optimization Toolbox] Optimization Toolboxのドキュメント

Optimization Toolboxのカテゴリが空だったのでとりあえずpdfドキュメントへのリンクを貼っておきます。

MATLAB:マニュアル:ドキュメンテーションセット:Optimization Toolbox:サイバネットシステム

一応WEBのドキュメントもありますが個人的にはキーワードで検索しづらいのがなぁ…と。
pdfの方がCtrl + FとかCtrl + Sで検索できて便利だと個人的には思うので。
あと印刷できるし。

それとこの辺り↓も参考になると思います。特に初心者には。
MATLAB:資料ライブラリ検索詳細:MATLABで行う最小化の基礎 - fminbnd関数とfminsearch関数 -:サイバネットシステム

で、まぁ結局読んでも上手くいかないと思うので,というか上手くいかなかったのであとは個々英語でなんとか打開してください…
The MathWorks - Support by Product - Optimization Toolbox

特にMATLAB CentralFile Exchangeに結構役立つものがあったりします。
The MathWorks - Support by Product - Optimization Toolbox - File Exchange
このエントリに付けられたタグ|MATLABOptimizationToolboxドキュメント