[2008/01/27]非線形方程式の解法なら、こちらのcsolveの記事を読むことをお勧めします。先日、Optimization Toolboxのfsolveを用いて，非線形方程式を解く方法を紹介しましたが，（
これと
これ），この方法は、Optimization Toolboxを購入している方のみしかできない方法です。
そこで，持ってない人が非線形方程式を解く方法の一例として，ニュートン・ラプソン法を実装してみました。
fsolveほど強力ではないですがそれなりに動くと思います…（急に書いてみたので保証はできませんが…）
ソースを書くと長くなるので肝になるところだけ書くと，こんな感じでヤコビ行列が指定された場合はそれを使用し，指定されない場合はここに書いた
jacobi_numeric_fast()で有限差分としてヤコビ行列を求め，ニュートン・ラプソン法を実行しています。
while exitflag==10
...
if strcmpi(opt.Jacobian, 'on')
[fs, jac] = func(x,varargin{:});
x = x - rat*jac¥fs;
else
fs = func(x,varargin{:});
jac = jacobi_numeric_fast(func, x, hh, varargin{:});
x = x - rat*jac¥fs;
end;
...
if fs < minfs minx = x; end;
% 終了条件のチェック
if max(abs(fs)) < opt.TolFun
exitflag = 1;
end;
...
end;
使い方は，一応fsolveに似せてみました。
こんな感じです。
x = newton_numeric(@fun,x0) … 初期値x0を用いて関数fun()で定義された非線形方程式を解き、解xを得る。
x = newton_numeric(@fun,x0,options) … 初期値x0を用い、オプションを構造体optionsで指定して関数fun()で定義された非線形方程式を解を解き、解xを得る。
x = newton_numeric(@fun,x0,options,fun()に渡すほかの引数...)
[x,fval,exitflag,output,jacobian] = newton_numeric(@fun,...)
x … 解
fval … 解xの時の最終的な関数funの値
exitflag … 計算アルゴリズムが停止した理由を示す整数（下記参照）
output … 最適化に関する情報を含む構造体（下記参照）
jacobian … 解xの時の関数funのヤコビ行列
基本的にfsolve()に使い方を似せたので，ヤコビ行列の指定など基本的な使い方については，
この記事と
この記事が参考になると思います。
以下，実行時に指定できるオプション，exitflag、outputの説明です。
まず，実行時に指定できるオプションについてです。
これはoptimset()でセットしたオプション構造体で問題ありません。
使用できるフィールドは以下の通り。
ここで，TolXとTolFunについてはMATLABのないでの一般的な扱いがわからなかったので，僕が独自に使用してます。
TolXとTolFunに関してはMATLAB標準の関数とは扱いが多少変わりますので，注意してください。| Display | 最小化あるいは計算する関数に関する情報の出力。（デフォルトは'final'） |
| MaxFunEvals | 許可する関数計算の最大回数 |
| MaxIter | 許可する繰り返しの最大回数 |
| Jacobian | 'on'の場合は、ユーザ定義のヤコビ行列を用い，'off'なら有限差分からヤコビ行列を近似する。（デフォルトは'off'） |
| TolFun | 関数値に関する終了許容範囲（maxFuncがTolFun以下なら解とみなす。） 注意！！ |
| TolX | xの変化に関する終了許容範囲 注意！！ |
次に出力として与えられる計算アルゴリズムが停止した理由を示す整数exitflagの説明を下記に書いておきます。
| 1 | 関数が解 xに収束したことを示します。 |
| 2 | 解xの変化が指定した許容値より小さくなったことを示します。 |
| 3 | 関数の変化が指定した許容値より小さくなった（あまり変化しなくなった）ことを示します。 |
| 0 | 関数計算の回数、あるいは繰り返し回数が最大回数に達していることを示します。 |
| -1 | 解が発散もしくは振動していることを示します。 |
また，最適化に関する情報を含む構造体outputの説明は下記です。
| iterations | 実行した繰り返し回数 |
| funcCount | 関数の評価回数 |
| message | 終了メッセージ |
| answer | 最終的な解x |
| Func | 目的関数の最終的な大きさ |
| maxFunc | 目的関数の要素のうち絶対が最大のものの大きさ |
| TolFun | 関数値に関する終了許容範囲（maxFuncがTolFun以下なら解とみなす。） |
| TolX | xの変化に関する終了許容範囲 |
以下に実際に例題を示します。
ここでは、fsolveで解いたとき（
これと
これ）と同様に下記の非線形方程式を解きます。関数も同じものを用いています。
まず，オプションを指定せず（ヤコビ行列も使用せず）解いた例です。
>> [x,fs,exitflag, output, jacobian] = newton_numeric(@funcsj, ones(3,1))
Iteration F(x) Diff. Pre-step Norm of Step
6 3.99456e-025 -1.8e-006 2.9e-013
x =
0.7331
1.8608
1.3641
fs =
1.0e-012 *
0.5498
-0.2633
-0.1669
exitflag =
1
output =
funcCount: 7
iterations: 6
TolX: 1.0000e-006
TolFun: 1.0000e-006
Func: 3.9946e-025
maxFunc: 5.4978e-013
answer: [0.7331 1.8608 1.3641]
message: '解が，許容範囲内に収束しました。'
jacobian =
1.4662 3.7216 0
1.8608 0.7331 -1.0000
1.3641 0 0.7331
次にオプションを指定し，ヤコビ行列も独自のものを用いた例です。
>> opt = optimset( 'Jacobian', 'on', 'Display', 'iter', 'TolX', 1e-12, 'TolFun', 1e-12);
>> [x,fs,exitflag, output, jacobian] = newton_numeric(@funcsj, ones(3,1), opt)
Iteration F(x) Diff. Pre-step Norm of Step
0 4 0 0
1 30 3.5 1.1
2 1.12809 -4.4 0.35
3 0.0157807 -0.94 0.056
4 8.02764e-006 -0.12 0.0013
5 3.13767e-012 -0.0028 8.3e-007
6 3.99456e-025 -1.8e-006 2.9e-013
x =
0.7331
1.8608
1.3641
fs =
1.0e-012 *
0.5498
-0.2633
-0.1669
exitflag =
1
output =
funcCount: 7
iterations: 6
TolX: 1.0000e-012
TolFun: 1.0000e-012
Func: 3.9946e-025
maxFunc: 5.4978e-013
answer: [0.7331 1.8608 1.3641]
message: '解が，許容範囲内に収束しました。'
jacobian =
1.4662 3.7216 0
1.8608 0.7331 -1.0000
1.3641 0 0.7331
このように簡単な例なら、先ほども書きましたがそれなりに動くと思います…
複雑な問題の場合にどうなるかはちょっと保証はできません。
それでもよろしいという方は下記のリンクから，zipファイルをダウンロードしてください。
newton_numeric()に同封で，jacobi_numeric_fast.mとデモのmファイルをいれてあります。
参考：Optimization Toolbox User's Guide - 最適化オプション参考：初心者によるMATLABメモ｜ 非線形方程式を解く(fsolve)参考：初心者によるMATLABメモ｜ 非線形方程式を解く2 (fsolve)参考：初心者によるMATLABメモ｜ ヤコビ行列（数値、Fast Code）ダウンロード：このエントリのmファイル
