初心者によるMATLABメモ

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

Google Adsenseのプライバシーポリシーについて明記したページへのリンクをサイドメニューに追加しました。

参考：初心者によるMATLABメモ | プライバシーポリシー

このブログも更新したいネタはいろいろあるんですが、仕事が忙しくなかなか更新できません・・・
機を見て一挙に更新したいと思います。

MATLABを使っていて、OSの環境変数を使用できると便利な場面が意外とあると思うのですが、そんな場合はgetenvを使用することで簡単に実現できます。

N = getenv('環境変数名')

例えば、自分はシミュレーションを複数の計算機でかけてるんですが、進捗状況メールがどの計算機から送られてきたのか混乱することも度々…
ということを繰り返した結果、今では進捗状況をメールするときに下記のようにサブジェクトにコンピュータ名を含めています。

subject = strcat(['[', getenv('COMPUTERNAME'), ']', 'テストメール']);

一応書いておくと、Windowsの場合は環境変数が見たければ、コマンドプロンプトで「set」コマンドを引数なしで実行すれば見れます。

参考：MATLAB 関数リファレンス - getenv

簡単な作業ログの保存などに使う関数diaryを拡張したdiary2()です。(mファイルのダウンロードはこちら)
基本的にはdiary()と機能、使い方は変わりません。

違う点は、diaryだと追加書き込みに自動的になってしまいますが、diary2だと新規書き込みを行うことができます。
また、開始時と停止時にタイムスタンプを自動的に書き込みます。
あとは現在のdiaryファイルの名前の表示、diaryのon/offの表示の程度の機能が追加されています。

使い方は基本的にはdiary()と同じです。

diary2　…　diary('on')と同じ
diary new　…　以後の出力をdiaryファイルに書き込みます。(ファイルが存在しても新規書き込み)
diary on　…　以後の出力をdiaryファイルに書き込みます。(ファイルが存在すれば追加書き込み)
diary off　…　書き込みを中止します。

diary stat　…　diaryモードのON/OFFを表示します。
diary fname　…　現在指定されているdiaryファイル名を表示します。

diary('ファイル名')　…　diaryファイル名を指定して書き込みを開始します。
diary('ファイル名','on')　…　diaryファイル名を指定し、追加書き込みを行うことができます。
diary('ファイル名','new')　…　diaryファイル名を指定し、新規書き込みを行うことができます。

※diary2('file')は、diary2('file','on')と同じです。

on/off/newあたりは普通のdiaryと使い方が同じなので、ここでは省略するとして、一応statとfnameだけ例を載せておきます。
一連の流れとして使用するとこんな感じです。

>> diary2('diary.txt','new')
>> diary2 fname
現在のdiaryのファイル名 ... [ diary.txt ]
>> diary2 stat
現在のdiaryの状態 ... [ on ]
>> diary2 off
>> diary2 stat
現在のdiaryの状態 ... [ off ]

最後にソースを載せておきます。
また、diary2()のヘルプ付きのmファイルのダウンロードはこちらから行えます。

function diary2(varargin) 

%% 引数のチェックと処理
error(nargchk(0,2,nargin));

opt = 0;
diary_path = get(0, 'DiaryFile');
mode = {'on','new','off','stat','fname'};

for ii = 1:nargin 
    in = varargin{ii};

    if ischar(in) 
        hit = find(strcmpi(lower(in), mode));
        if isempty(hit) 
            if ~strcmpi(diary_path,in) 
                diary2 off;
                diary_path = in;
            end;
        else 
            opt = hit(1);
        end;
    else 
        error('引数は文字列である必要があります');
    end;
end;

if opt==0 
    if strcmpi(get(0,'Diary'),'on') 
        opt=3;
    else 
        opt=1;
    end;
end;
%% 


%% 実際の処理
switch opt 
    % 書き込み開始 
    case {1,2} 

        % diaryファイルへの書き込み開始メッセージ書き込み 
        if strcmpi(get(0,'Diary'),'off') 
            if opt==1 
                perm = 'a+';
            else 
                perm = 'w+';
            end;
            fid = fopen(diary_path,perm);
            fprintf(fid, '%s¥n', '%%% =================  diary 書き込み開始 ==================== ');
            fprintf(fid, '%s%s¥n¥n', '+ 現在時間 : ', datestr(now,31));
            fclose(fid);
        end;

        % diary on
        diary(diary_path);
    % 書き込み停止 
    case 3
        % diaryファイルへの書き込み終了メッセージ書き込み 
        if strcmpi(get(0,'Diary'),'on') 
            diary off;

            fid = fopen(diary_path,'r+');

            body = {};
            frewind(fid);
            while 1 
                tline = fgetl(fid);
                if ~ischar(tline), break, end;
                body{length(body)+1} = tline;
            end;

            body{end} = [];

            frewind(fid);
            for ii=1:length(body)
                fprintf(fid, '%s¥n', body{ii});
            end;
            fprintf(fid, '%s%s¥n', '+ 現在時間 : ', datestr(now,31));
            fprintf(fid, '%s¥n', '%% ========== / diary 書き込み停止 ==========');

            fclose(fid);
        end;
    % ステータス 
    case 4 
        disp(sprintf('現在のdiaryの状態 ... [ %s ]', get(0,'Diary')));
    % ファイル名表示 
    case 5 
        disp(sprintf('現在のdiaryのファイル名 ... [ %s ]', get(0,'DiaryFile')));
    otherwise
        help diary2;
end

参考：MATLAB Function Reference - diary
ダウンロード：このエントリのmファイル　左クリックでお願いします。　（リファラを送信しない場合は，アクセス拒否になってしまいます）

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

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からメールを送信することができます。
これが意外と便利なので紹介しておきます。
例えば、どれくらいかかるかわからない長い計算をかけている場合、「終わったかな？」とわざわざチェックするのは面倒です。
そこで、僕は計算が終わったらメールが来るようにしています。

以下、MATLABからメールを送信する方法です。
まず、基本的なSMTPサーバ（送信用のメールサーバ）の設定を行います。
ここらへんはstartup.mに書いておくと便利でしょう。ただし、パスワードも記入する場合には取り扱いに注意してください。

% SMTPサーバの設定
setpref( 'Internet', 'SMTP_Server', 'SMTPサーバ');
% SMTPサーバでのユーザ名（必要なら）
setpref( 'Internet', 'SMTP_Username', 'ユーザ名');
% SMTPサーバでのパスワード（必要なら）
setpref( 'Internet', 'SMTP_Password', 'パスワード');
% 自分のアドレス（Fromになる）
setpref('Internet', 'E_mail', '自分のアドレス');

あとは下記のようにsendmail関数を使用すればメールを送ることが出来ます。

sendmail('相手のメールアドレス','題名','メッセージ')

また、添付ファイルをつけたい場合は以下です。

sendmail('相手のメールアドレス','題名','メッセージ','添付ファイルのパス')

その他、複数の人に送る場合、複数の添付ファイルを付けたい場合などはこちらを参考にしてください。

参考：MATLAB：マニュアル:サイバネットシステム - MATLAB Function Reference - sendmail
参考：MATLAB：マニュアル:サイバネットシステム - MATLAB Programming - E メールの送信
参考：初心者によるMATLABメモ｜ 起動時にデフォルトカラーなどを設定する　（起動時にスクリプトを実行）