FC2ブログ

初心者によるMATLABメモ

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

 

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

[グラフ、プロット] プロット(グラフ)で目盛りのラベルを指定する

MATLABでプロット(グラフ)で目盛りのラベルを手動で設定する方法です。
目盛りの刻みを手動で設定するときと同様に、特定の関数は用意されておらず,直接Axesオブジェクトのプロパティをset()で弄る必要があります。
目盛りのラベルを手動で設定
set(gca,'XTickLabel', 目盛りのラベルとして使う文字行列);
set(gca,'YTickLabel', 目盛りのラベルとして使う文字行列);
set(gca,'ZTickLabel', 目盛りのラベルとして使う文字行列);

※目盛りのラベルとして使う文字行列としては下記のフォーマットが使用できます。
{'AAA';'BBB';'CCC'} … 文字列を収めたセル配列
['A';'B';'C'] … 文字列配列
'A|B|C|D|E|F|G' … 「|」で区切った文字列
[100,-100,300] … 数値配列

まず、棒グラフのX軸に単にラベルを指定した(駄目な)例です。
実行結果を見てもらうとわかると思うのですが、この場合ラベルが意図しない位置についてしまっています。
data = [1,2,5,10,4,6];
figure('InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [200 100 320 160]);
hold all;
bar(data);
set(gca, 'XTickLabel',{'1月';'2月';'3月';'4月';'5月';'6月'});



このようになるのを避けるには下の例ようにラベルを指定すると同時に目盛りの刻みも設定してやる必要があります。
tcksX={'1月';'2月';'3月';'4月';'5月';'6月'};
set(gca, 'XTickLabel',tcksX ,'XTick',1:length(tcksX));



さらにY軸に指定のラベルを等間隔で指定する方法です。
一度先にylim()でy軸の範囲を取得してから,それをlinspace()で当分割してY軸の目盛りを指定しています。
tcksX={'1月';'2月';'3月';'4月';'5月';'6月'};
tcksY={'A';'B';'C';'D';'E'};
yrange = ylim();
set(gca, 'XTickLabel',tcksX ,'XTick',1:length(tcksX), 'YTickLabel',tcksY , 'YTick', linspace(yrange(1), yrange(2), length(tcksY)));



最後にこれは駄目な例ですが、このように目盛りの数と合っていないラベルを指定すると、ラベルがループするようになっています。
set(gca, 'XTickLabel','1月|2個目');



参考:MATLAB Function Reference - Axes Properties
参考:MATLAB Function Reference - get
参考:MATLAB Function Reference - set
スポンサーサイト



このエントリに付けられたタグ|MATLABグラフプロット目盛り目盛りのラベル目盛りのテキストAxes

[Graphics] 現在のカラーマップを関数として保存する(save_colormap)

MATLABでFigureカラーマップエディタでカラーマップを弄ったときにそのカラーマップをとっておきたいというときがあると思います。
そんなときに使える現在のカラーマップを独自の関数とし保存する関数save_colormapを紹介します。

コードは下記になります。
やっていることは現在のカラーマップを取得して,テキストファイルに書き出しているだけです。
このコードをコピペするより,エラーチェックなどをもう少ししているので下のリンクからmファイルをダウンロードした方が良いかと思います。(このエントリのmファイルへのリンク
function save_colormap(fname, dpath) 

if nargin<1 || isempty(fname) fname = 'mycolormap'; end;
if nargin<2 dpath = pwd; end;

if ~ischar(dpath) error('第2引数は文字列である必要があります。'); end; 

% ディレクトリのパスのチェック
[dirname, filename, ext, versn] = fileparts(dpath);
[num,pathinfo] = fileattrib(fullfile(dirname, filename));
if pathinfo.directory~=1 error('第2引数がディレクトリのパスではありません。'); end;
fpath = fullfile(pathinfo.Name, strcat([fname, '.m']));

% 現在のカラーマップの取得
fgh = get(0,'CurrentFigure');
if isempty(fgh) error('カレントのFigureがありません。'); end;
ccl = get(fgh,'colormap');

% カラーマップの書き込み
fh = fopen(fpath, 'w');
fprintf(fh, '%s%s%s\n%s', 'function CMAP = ', fname, '()', 'cols = [' );
fprintf(fh, '%d ',ccl);
fprintf(fh, '%s\n%s\n%s', '];', 'CMAP = reshape(cols, floor(length(cols)/3), 3);');
fclose(fh);

一応簡単な使い方です。
関数名を省略すると関数名「mycolormap」として,ディレクトリを省略すると現在のディレクトリに保存します。
save_colormap('関数名', '保存するディレクトリ')

以下、使用例です。
まず,曲面を表示します。
>> [x,y] = meshgrid([-2:.2:2]);
Z = x.*exp(-x.^2-y.^2);
surf(x,y,Z,gradient(Z))
colorbar;



続いて、Figureウインドウの編集(Edit)メニューからカラーマップ(Colormap)を選択し、カラーマップエディタを起動してカラーマップを変更します。




現在のカラーマップをmycolormap.mとして保存します。
>> save_colormap()

カラーマップをspringに変更します。
>> colormap spring;



カラーマップを先ほど保存したカラーマップに戻してみます。
>> colormap mycolormap;



とりあえず僕が書いたmファイルは下記リンクに置いてあります。

参考:MATLAB Function Reference - exist
参考:MATLAB Function Reference - colormap
参考:MATLAB Function Reference - colormapeditor
このエントリのmファイル(save_colormap.m)
このエントリに付けられたタグ|MATLABグラフプロットカラーマップ保存save_colormapcolormap

[グラフ、プロット] MATLABからエクセル(Microsoft Excel)にグラフを描く

たまたま他のものをFile Exchangeで探していたら、matlabからExcelにグラフを描くxlschartってのを見つけました。
僕個人としてはMATLABでグラフを描いた方が綺麗だし、融通が利くのでいいと思うんですが、あまりMATLABを使わないのにプロットのカスタマイズまで覚えるのは…という人も多そうで需要がありそうなので載せておきます。

簡単な使い方と引数の説明はこんな感じです。詳しくはヘルプを参照してください。(日本語のヘルプはこちらの記事をご覧ください。)
xlschart(column_labels,data,chartype,chart_title)
xlschart(column_labels,data,chartype,chart_title,filename)
xlschart(column_labels,data,chartype,chart_title,xtitle,ytitle)
xlschart(column_labels,data,chartype,chart_title,filename,sheetname)
xlschart(column_labels,data,chartype,chart_title,xtitle,ytitle,filename)
xlschart(column_labels,data,chartype,chart_title,xtitle,ytitle,filename,sheetname)

※引数について
column_labels : 各列のラベル(文字列のセル配列で指定)。
data : データ(数値行列)。
chartype : グラフの種類の整数値。または,文字列で指定(例えば'3DLine')。
chart_title : グラフのタイトルの文字列とシート名.。
filename : エクセルファイルの名前.(これが指定された場合は自動的に保存され,Excelは立ち上がらない。)。
sheetname:書き込むエクセルのシート名。
xtitle : X軸のラベル(column_labelsに指定したセル行列に含まれている必要がある)。
ytitle : Y軸のラベル(column_labelsに指定したセル行列に含まれている必要がある)。

以下、実際に使った例を2つほど。
まずは棒グラフを描いた例です。
画像を見てもらえば分かるとおり,Excel内にグラフが描かれます。(注:色、フォントサイズなどはExcelで弄ってあります。)
columns = {'データX','データY'};
data = [0:6;rand(1,7)]';
xlschart(columns, data, 1, 'MATLABからExcelにグラフを描く - 棒グラフ', 'データX', 'データY')



次は、株価チャートをを描いてみた例です。(注:色、フォントサイズなどはExcelで弄ってあります。)
columns = {'日付', '出来高', '始値', '高値', '安値', '終値'};
nn=7;

days=datestr(now+[1:nn],'yyyymmdd');
days=str2num(days);
vv=floor(50000*rand(nn,1));
hh=48+rand(nn,1);
ll=44+rand(nn,1);
oo=46+2*(0.5-rand(nn,1));
cc=46+2*(0.5-rand(nn,1));

data = [days,vv,oo,hh,ll,cc]
xlschart(columns, data, 52, 'MATLABからExcelにグラフを描く - 株価チャートの例');




参考:初心者によるMATLABメモ| エクセル(Microsoft Excel)でのグラフに対応するMATLABの関数
参考:MATLAB Central File Exchange - xlschart
参考:初心者によるMATLABメモ| MATLABからExcelにグラフを描くxlschartの日本語ヘルプ
このエントリに付けられたタグ|MATLABグラフExcelエクセルxlschart

[グラフ、プロット] MATLABからExcelにグラフを描くxlschartの日本語ヘルプ

この記事で紹介したMATLABからエクセル(Microsoft Excel)にグラフを描くxlschartの日本語ヘルプを書いてみました。

使い方としては,ダウンロードしたxlschart.mを開き、下記の2行目から59行目あたりの下に,僕が書いた日本語ヘルプを付け加えればよいかと思います。
結構長いのでコピペしづらいかなと思ったので一応ヘルプだけを書いたテキストファイルを用意しました。必要な方はこちらよりダウンロードしてください。
 1:function xlschart(titles,m,chartype,chart_title,varargin)

 2:%XLSCHART Creates graphs in Microsoft Excel.



59:% See also XLSREAD, XLSFINFO, XLSWRITE, XLSCELL, XLSHEETS, , CPTXT2XLS, MSOPEN

以下、僕が書いたxlschartの日本語ヘルプです。
% XLSCHART    Microsoft Excelの中にグラフを描く.

% xlschart(column_labels,data,chartype,chart_title)
% xlschart(column_labels,data,chartype,chart_title,filename)
% xlschart(column_labels,data,chartype,chart_title,xtitle,ytitle)
% xlschart(column_labels,data,chartype,chart_title,filename,sheetname)
% xlschart(column_labels,data,chartype,chart_title,xtitle,ytitle,filename)
% xlschart(column_labels,data,chartype,chart_title,xtitle,ytitle,filename,sheetname)
%
% xlschart  :Microsoft Excelにデータのカラムヘッダと行列形式でデータを書き込み,
%            さらに2つ以上の行のデータを用いてグラフを描きます。
%
% 引数について
%       column_labels:  各列のラベル(文字列のセル配列で指定)。
%       data:           データ(数値行列)。
%       chartype:       下のリストに表示されたグラフの種類の整数値。
%                       または,文字列で指定(例えば'3DLine')。
%       chart_title:    グラフのタイトルの文字列とシート名.。
%       filename:       エクセルファイルの名前.(これが指定された場合は自動的に保存され,Excelは立ち上がらない。)。
%       sheetname:      書き込むエクセルのシート名。
%       xtitle:         X軸のラベル(column_labelsに指定したセル行列に含まれている必要がある)。
%       ytitle:         Y軸のラベル(column_labelsに指定したセル行列に含まれている必要がある)。

% CHARTYPEに指定できる整数値と名称とグラフの種類:
%    1- ColumnClustered  -  集合縦棒グラフ
%    2- ColumnStacked  -  積み上げ縦棒グラフ
%    3- ColumnStacked100  -  100%積み上げ縦棒グラフ
%    4- 3DColumnClustered  -  3Dの集合縦棒グラフ
%    5- 3DColumnStacked  -  3Dの積み上げ縦棒グラフ
%    6- 3DColumnStacked100  -  3Dの100%積み上げ縦棒グラフ
%    7-3DColumn  -  3Dの縦棒グラフ(系列間と項目間の比較)
%    8-BarClustered  -  集合横棒グラフ
%    9-BarStacked  -  積み上げ横棒グラフ
%    10-BarStacked100  -  100%積み上げ横棒グラフ
%    11-3DBarClustered  -  3Dの集合横棒グラフ
%    12-3DBarStacked  -  3Dの積み上げ横棒グラフ
%    13-3DBarStacked100  -  3Dの100%積み上げ横棒グラフ
%    14-Line  -  折れ線グラフ
%    15-LineStacked  -  積み上げ折れ線グラフ
%    16-LineStacked100  -  100%積み上げ折れ線グラフ
%    17-LineMarkers  -  マーカー付きの折れ線グラフ
%    18-LineMarkersStacked  -  マーカー付きの積み上げ折れ線グラフ
%    19-LineMarkersStacked100  -  マーカー付きの100%積み上げ折れ線グラフ
%    20-3DLine  -  3D折れ線グラフ
%    21-Pie  -  円グラフ
%    22-3DPie  -  3D円グラフ
%    23-PieOfPie  -  補助円グラフ付き円グラフ
%    24-PieExploded  -  分割円グラフ
%    25-3DPieExploded  -  3Dの分割円グラフ
%    26-BarOfPie  -  補助縦棒付き円グラフ
%    27-XYScatter  -  散布図
%    28-XYScatterSmooth  -  平滑線付きの散布図(マーカー有)
%    29-XYScatterSmoothNoMarkers  -  平滑線付きの散布図(マーカー無)
%    30-XYScatterLines  -  折れ線付きの散布図(マーカー有)
%    31-XYScatterLinesNoMarkers  -  折れ線付きの散布図(マーカー無)
%    32-Area  -  面グラフ
%    33-AreaStacked  -  積み上げ面グラフ
%    34-AreaStacked100  -  100%積み上げ面グラフ
%    35-3DArea  -  3Dの面グラフ
%    36-3DAreaStacked  -  3Dの積み上げ面グラフ
%    37-3DAreaStacked100  -  3Dの100%積み上げ面グラフ
%    38-Doughnut  -  ドーナツグラフ
%    39-DoughnutExploded  -  分割ドーナツグラフ
%    40-Radar  -  レーダーチャート
%    41-RadarMarkers  -  マーカー付きのレーダーチャート
%    42-RadarFilled  -  塗りつぶしレーダーチャート
%    43-Surface  -  等高線グラフ
%    44-SurfaceWireframe  -  ワイヤフレーム等高線グラフ
%    45-SurfaceTopView  -  上から見た等高線グラフ
%    46-SurfaceTopViewWireframe  -  上から見たワイヤフレーム等高線グラフ
%    47-Bubble  -  バルブチャート
%    48-Bubble3DEffect  -  3Dバルブチャート
%    49-StockHLC  -  株価チャート(高値 - 安値 - 終値)
%    50-StockOHLC  -  株価チャート(始値 - 高値 - 安値 - 終値)
%    51-StockVHLC  -  株価チャート(出来高 - 高値 - 安値 - 終値)
%    52-StockVOHLC  -  株価チャート(出来高 - 始値 - 高値 - 安値 - 終値)
%
% 例:
%      column_labels = {'1st','2nd','3rd','4th','5th','6th','7th','8th','9th','10th'};
%      m = magic(10);
%      xlschart(column_labels,data,'LineStacked100','My Title')
%      xlschart(column_labels,data,17,'My Title','sample.xls','Sheet2');
%      xlschart(column_labels,data,'Line','My Title','1st','2nd','sample.xls','Sheet2');
%      xlschart(column_labels,data,20,'My Title','1st','8th','sample.xls','Sheet2');
%      xlschart(column_labels,data,'XYScatterSmooth','My Title','10th','9th','sample.xls','Sheet2');
%      xlschart(column_labels,data,30,'My Title','2nd',{'7th','6th'},'sample.xls','Sheet2');
%      xlschart(column_labels,data,29,'My Title',{'2nd','5th'},{'7th','6th'},'sample.xls','Sheet2');
%
%   参考 XLSREAD, XLSFINFO, XLSWRITE, XLSCELL, XLSHEETS, , CPTXT2XLS, MSOPEN

参考:MATLAB Central File Exchange - xlschart
参考:初心者によるMATLABメモ| MATLABからエクセル(Microsoft Excel)にグラフを描く
このエントリのxlschartの日本語ヘルプのファイル
このエントリに付けられたタグ|MATLABグラフExcelエクセルxlschart日本語ヘルプ

[基本] 2進数、16進数と10進数の相互変換

MATLABで使える2進数、16進数と10進数を相互に変換する関数です。
bin2dec() … 2進数(文字列)を10進数に変換。
dec2bin() … 10進数整数を2進数(文字列)に変換。
hex2dec() … 16進数(文字列)を10進数に変換。
dec2hex() … 10進数整数を16進数(文字列)に変換。
dec2hex(bin2dec()) … 2進数(文字列)を16進数(文字列)に変換。
dec2bin(hex2dec()) … 16進数(文字列)を2進数(文字列)に変換。

hex2num … 16進数(文字列)を倍精度数値(double)に変換。

以下、簡単に使用例。

2進数(文字列)→10進数
>> bin2dec('11111111')
ans =
   255

10進数→2進数(文字列)
>> dec2bin(255)
ans =
11111111
>> ischar(dec2bin(255))    % 文字列か確認
ans =
     1

16進数(文字列)→10進数
>> hex2dec('FF')
ans =
   255

10進数→16進数(文字列)
>> dec2hex(255)
ans =
FF
>> ischar(dec2hex(255))    % 文字列か確認
ans =
     1

2進数(文字列)→16進数(文字列)
2進数を16進数に変換する場合には,一度10進数に変換してやる必要があります。(逆も同様)
>> dec2hex(bin2dec('11111111'))
ans =
FF

16進数(文字列)→2進数(文字列)
>> dec2bin(hex2dec('FF'))
ans =
11111111

この辺りの変換をよく使う人は,こんな関数hex2bin()とbin2hex()を定義しておいた方がいいかもしれません。
function h = hex2bin(b)
h = dec2hex(bin2dec(b));

function b = bin2hex(h)
b = dec2bin(hex2dec(h));


参考:MATLAB Function Reference データタイプ データタイプの変換

このエントリに付けられたタグ|MATLAB2進数16進数10進数変換

[Graphics] MATLABの既定のカラーが省略名からRGB値のカラーベクトルに変換する

MATLABでは,色については通常はカラーベクトルで指定しますが,ここに書いてあるように,8個の既定のカラーについては青なら'b'といったような省略名または完全名での指定もすることが出来ます。
RGB値カラーベクトル省略名完全名
[1 1 0]yyellow
[1 0 1]mmagenta 
[0 1 1]ccyan
[1 0 0]rred
[0 1 0]ggreen
[0 0 1]bblue
[1 1 1]wwhite
[0 0 0]kblack

ここに載せている関数color_shortcut()は,これらのMATLABの8個の既定のカラーの省略名または完全名をRGB値のカラーベクトルに変換する関数です。
この関数は,実用な用途は低いかもしれませんが,自作の関数のサブ関数として使えると思いますので、参考までに紹介しておきます。
一応例によってヘルプ付きのmファイルはこちらからダウンロードできます。
function COL = color_shortcut(v) 

colname = {'b','y','m','c','r','g','w','k', ...
    'blue','yellow','magenta','cyan','red','green','white','black', ...
    };

boo = (ischar(v) && any(strcmpi(v,colname))) || ...
    (isrealvec(v) && numel(v)==3 && all(v>=0) && all(v<=1));
if boo~=1 error('色の指定が不正です'); end;

COL = zeros(1,3);
if ischar(v) 
    switch mod(find(strcmpi(lower(v), colname)),8) 
        case 1    % 'b' or 'blue'
            COL = [0 0 1];
        case 2    % 'y' or 'yellow'
            COL = [1 1 0];
        case 3    % 'm' or 'magenta'
            COL = [1 0 1];
        case 4    % 'c' or 'cyan'
            COL = [0 1 1];
        case 5    % 'r' or 'red'
            COL = [1 0 0];
        case 6    % 'g' or 'green'
            COL = [0 1 0];
        case 7    % 'w' or 'white'
            COL = [1 1 1];
        case 0    % 'k' or 'black'
            COL = [0 0 0];
    end;
else 
    COL = v;
end

一応実行例。
この例では,青色のカラーベクトル[0 0 1]をCとして取得しています。
>> C = color_shortcut('b')
C =
     0     0     1


参考:MATLAB Function Reference - ColorSpec
このエントリのmファイル(color_shortcut())
このエントリに付けられたタグ|MATLABGraphics規定のカラー省略名完全名color_shortcut

[Geometry] 単位ベクトル2 (unitvector)

前に書いた方法だと列ベクトルとしてひとつのベクトルを表した場合にしか対応していなかったのを拡張してみました。
以下,簡単な使い方,ソース,使用例の順に示します。
例によってヘルプ付きmファイルが欲しい方は,こちらからダウンロードできます。
ダウンロードファイルはテキストファイルになっているのでunitvector.m.txtからunitvector.mに名前を変更してください。

まずは簡単な使い方。
unitvector(V) … 列ベクトルとしてひとつのベクトルが表されたベクトル配列Vの各列の単位ベクトルを返します。
unitvector(V, DIM)は単位化する次元を指定したベクトル配列Vの各列の単位ベクトルを返します。

※unitvector(V)は、unitvector(V, 1)と同じです。

続いて,ソースです。
多少のエラー処理は一応しましたが、十分でないです。
function UV = unitvector(V,dim)

error(nargchk(1, 2, nargin));

if nargin==1 dim=1; end;

if ~isreal(V) 
    error('ベクトルは実数配列である必要があります。'); 
end;

if ~isequal(dim,floor(dim)) || numel(dim)~=1 || ndims(V)    error('次元の引数は、インデックス付け範囲の中で、正の整数のスカラでなければなりません。');
end;

sz = size(V);
sz2 = ones(1,ndims(V));
sz2(dim) = sz(dim);

vnorm = sqrt(sum(V.^2,dim));
vnorm = repmat(vnorm,sz2);
vnorm = reshape(vnorm,sz);

if any(vnorm==0) 
    warning('ゼロ割りです。NaNもしくはInfになっているベクトルがあります。');
end;

UV = V./vnorm;

最後に使用例をいくつか示します。
まずは1つのベクトルを列ベクトルで表したベクトル集合をそれぞれ単位ベクトルにします。
>> A=[[1;1;0],[1;0;0]]
A =
     1     1
     1     0
     0     0
>> unitvector(A)
ans =
    0.7071    1.0000
    0.7071         0
         0         0

続いて、1つのベクトルを行ベクトルで表したベクトル集合をそれぞれ単位ベクトルにします。
>> A = [[1,1,0];[1,0,0]]
A =
     1     1     0
     1     0     0
>> unitvector(A,2)
ans =
    0.7071    0.7071         0
    1.0000         0         0

次はxyzをそれぞれ次元を変えて表現した場合です。
>> A = reshape([[1,1,0];[1,0,0]],1,2,3)
A(:,:,1) =
     1     1
A(:,:,2) =
     1     0
A(:,:,3) =
     0     0
>> unitvector(A,3)
ans(:,:,1) =
    0.7071    1.0000
ans(:,:,2) =
    0.7071         0
ans(:,:,3) =
     0     0


参考:初心者によるMATLABメモ| 単位ベクトル
ダウンロード:このエントリのmファイル
このエントリに付けられたタグ|MATLABGeometry幾何単位ベクトルunitvector修正版

[Graphics] カラーベクトルから手動でグラデーションを作る

カラーベクトル間を線形補間し,指定した色数のグラデーションを作成する関数です。
やっていることは書いたそのままなので説明はしませんが,簡単な使い方,使用例の順に示します。(ダウンロードはこちらから)
ちなみにこの関数は,ここにあるlinspacem()を内部で使用しているので,linspacem()も合わせてダウンロードし,同じフォルダ(ディレクトリ)またはパスの通っているところに置かないとエラーになります。

まずは、簡単な使用方法です。
ez_gradation(CVEC) … 複数のカラーベクトルCVECからグラデーションを作成し,色数64のカラーマップを作成する。
ez_gradation(CVEC, CNUM) … 複数のカラーベクトルCVECからグラデーションを作成し,色数CNUMのカラーマップを作成する。

※MATLABの仕様により、カラーマップの長さは任意(MS-WindowsとMacintosh上では256まで)です。

続いて簡単な使用例。
まずはこんな感じの曲面を生成。
[X,Y] = meshgrid(0:30);
Z = exp(X/10);
surf(X,Y,Z);
axis equal;
colorbar;
set(gcf, 'InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [200 100 380 260]);



緑[0,1,0]と赤[1,0,0]からなるカラーベクトルからグラデーションを作成し,適用してみます。
C=ez_gradation([[0,1,0];[1,0,0]]);
colormap(C);



色数を4に減色した緑[0,1,0]と赤[1,0,0]からなるカラーベクトルからグラデーションを作成し,適用した例です。
C=ez_gradation([[0,1,0];[1,0,0]], 4);
colormap(C);



最後に乱数を利用して作成した色数128のカラーマップを適用してみた例です。
C=ez_gradation(sin(pi*rand(5,3)/2), 128);
colormap(C);



別にカラーマップエディタで作ってもいいんですが,いちいち作るのも面倒くさい場面もあると思うので紹介してみました。

参考:初心者によるMATLABメモ| 等間隔に並んだ点列を取得する(linspacem)
ダウンロード:このエントリのmファイル(ez_gradation.m) ※この関数の動作は,linspacem()が必要です。
このエントリに付けられたタグ|MATLABGeometryカラーマップcolormapグラデーションlinspacemez_gradation

[数値解析] 双線形補間(バイリニア補間)のためのウェイト(重み)

双線形補間(バイリニア補間)のためのウェイト(重み)を生成する関数bilinear_weight()です。
簡単な使い方,コード、使用例の順に示します。(ダウンロードはこちらから)

まずは、簡単な使用方法です。
W = bilinear_weight(N) … 双線形補間のためのNxNマトリックスで表されるウェイトの長さ4のセル配列Wを生成します。
W = bilinear_weight(N,M) … 双線形補間のためのNxMマトリックスで表されるウェイトの長さ4のセル配列Wを生成します。
[W1,W2,W3,W4] = bilinear_weight(N,M) … 双線形補間のためのNxMマトリックスで表されるウェイトをそれぞれW1,W2,W3,W4として生成します。

この関数bilinear_weight()が出力しているウェイトを用いて,補間対象が下記の形で表される場合に,a*w1+b*w2+c*w3+d*w4またはa*w{1}+b*w{2}+c*w{3}+d*w{4}とすることで双線形補間を実行することができます。
|a  c|
|b  d|

実際にコードはこんな感じです。
function [w1,w2,w3,w4] = bilinear_weight(yn, xn) 

if nargin == 0 yn = 100; end;
if nargin == 1 xn = yn; end;

if ~isequal([yn,xn],floor([yn,xn])) || numel(xn)~=1  || numel(yn)~=1  
  error('引数は正のスカラである必要があります。');
end;

xn = floor(xn);
yn = floor(yn);

w1 = [linspace(1,0,yn)']*linspace(1,0,xn);
w2 = w1(end:-1:1,:);
w3 = w1(:,end:-1:1);
w4 = w1(end:-1:1,end:-1:1);

if any(nargout==[0 1]) 
  w1 = {w1,w2,w3,w4};
elseif nargout~=4 
  error('出力は1または4だと想定されています。');
end;

つづいて使用例です。
双線形補間のための3x4マトリックスで表されるウェイトをそれぞれw1,w2,w3,w4として生成した例。
>> [w1,w2,w3,w4]=bilinear_weight(3,4)
w1 =
    1.0000    0.6667    0.3333         0
    0.5000    0.3333    0.1667         0
         0         0         0         0

w2 =
         0         0         0         0
    0.5000    0.3333    0.1667         0
    1.0000    0.6667    0.3333         0

w3 =
         0    0.3333    0.6667    1.0000
         0    0.1667    0.3333    0.5000
         0         0         0         0

w4 =
         0         0         0         0
         0    0.1667    0.3333    0.5000
         0    0.3333    0.6667    1.0000

さらに、実際にこのウェイト[w1,w2,w3,w4]を用いて簡単な双線形補間を実行した例を下記に示します。
>> A=[[1;2],[3;4]]
A =
     1     3
     2     4

>> A(1)*w1+A(2)*w2+A(3)*w3+A(4)*w4
ans =
    1.0000    1.6667    2.3333    3.0000
    1.5000    2.1667    2.8333    3.5000
    2.0000    2.6667    3.3333    4.0000

最後にzを乱数で指定した入力メッシュを作成して、バイリニア補間した曲面を生成してみます。
まず入力メッシュです。
nn=3;
[px,py] = meshgrid(0:nn);
pz=zeros(size(px));
pz(2:end-1,2:end-1)=2*(0.5-rand([size(px)-2]));
pp = zeros(nn+1,nn+1,3);
pp(:,:,1) = px;
pp(:,:,2) = py;
pp(:,:,3) = pz;
psz = size(pp);

figure('InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [200 100 380 320]);
hold on;
mesh(pp(:,:,1),pp(:,:,2),pp(:,:,3), 'Marker', 'o', 'FaceAlpha', 0, 'EdgeColor', [1 0 0]);
surf(pp(:,:,1),pp(:,:,2),pp(:,:,3));
axis equal;
view(43,23)
hold off;


続いて、この入力メッシュをバイリニア補間した曲面を生成し,表示します。
上下の画像からちゃんと双線形補間されていることが確認できると思います。
sdiv = 10;
qsz = floor([sdiv*psz(1),sdiv*psz(2),psz(3)]);
qq = zeros(qsz);

ph = floor(linspace(1,qsz(1),psz(1)));
sdivh = diff(ph) + 1;

pw = floor(linspace(1,qsz(2),psz(2)));
sdivw = diff(pw) + 1;

% バイリニア補間
for jj = 1:(psz(1)-1) 
  for ii=1:(psz(2)-1) 
    w = bilinear_weight(sdivh(jj),sdivw(ii));
    for kk=1:psz(3) 
      qq(ph(jj):ph(jj+1) , pw(ii):pw(ii+1),kk) = pp(jj,ii,kk)*w{1} + pp(jj+1,ii,kk)*w{2} + pp(jj,ii+1,kk)*w{3} + pp(jj+1,ii+1,kk)*w{4};
    end;
  end;
end;

figure('InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [200 100 380 320]);
hold on;
mesh(pp(:,:,1),pp(:,:,2),pp(:,:,3), 'Marker', 'o', 'FaceAlpha', 0, 'EdgeColor', [1 0 0]);
surf(qq(:,:,1),qq(:,:,2),qq(:,:,3));
axis equal;
view(43,23)
hold off;



また、この関数bilinear_weightを用いて画像を拡大・縮小する例をこちらに掲載してます。こちらも参考にしてください。

参考:初心者によるMATLABメモ| 双線形補間(バイリニア補間)を用いた画像の拡大・縮小
ダウンロード:このエントリのmファイル(bilinear_weight.m)
このエントリに付けられたタグ|MATLAB数値解析線形補間双線形補間バイリニア補間bilinear_weight重みウェイト

[Graphics] 双線形補間(バイリニア補間)を用いた画像の拡大・縮小

こちらの記事に載せた双線形補間(バイリニア補間)のためのウェイト(重み)を生成する関数bilinear_weight()を利用して画像を拡大・縮小してみる例を載せてみたいと思います。(bilinear_weight()をダウンロードする場合はこちらから)

実際には,画像をバイリニア補間により拡大縮小する関数img_zoon_bilinear()を下記のように作成して、JPEGを拡大・縮小してみました。
この関数はあくまでもこの実験用に作ったもので,エラーチェックなど一切していないので気をつけてください。
function img_zoon_bilinear(i_filename, o_filename, ratio) 

% jpeg読み込み
i_img = imread(i_filename);
i_img = double(i_img)/255;

i_sz = size(i_img);
if numel(i_sz)<3 i_sz(3)=1; end;

% 出力画像の準備
o_sz = floor([ratio*i_sz(1),ratio*i_sz(2),i_sz(3)]);
o_img = zeros(o_sz);

oh = floor(linspace(1,o_sz(1),i_sz(1)));
sdivh = diff(oh) + 1;

ow = floor(linspace(1,o_sz(2),i_sz(2)));
sdivw = diff(ow) + 1;

wdef = bilinear_weight(round(ratio)+1,round(ratio)+1);

% バイリニアで補間
for jj = 1:(i_sz(1)-1) 
    for ii=1:(i_sz(2)-1) 
        if all([sdivh(jj),sdivw(ii)] == round(ratio)+1) 
            w = wdef;
        else 
            w = bilinear_weight(sdivh(jj),sdivw(ii));
        end;

        for kk=1:i_sz(3) 
            o_img(oh(jj):oh(jj+1) , ow(ii):ow(ii+1),kk) = i_img(jj,ii,kk)*w{1} + i_img(jj+1,ii,kk)*w{2} + i_img(jj,ii+1,kk)*w{3} + i_img(jj+1,ii+1,kk)*w{4};
        end;
    end;
end;

% 書き出し
imwrite(o_img, o_filename);

以下、実行例です。
まず、カラー画像に対して実行してみた例です。
これ↓が元画像です


以下が、倍率を変えて拡大・縮小してみた例です。
(50%) (80%) (120%) (200%)

また、200%のものをさらに50%にして元の大きさにしたものも載せておきます。


また、グレースケールのモノクロ画像に対しても、同様に実行することができます。
 元画像

(50%) 
(200%)


参考:MATLAB Function Reference - image
参考:MATLAB Function Reference - imread
参考:MATLAB Function Reference - imwrite
参考:初心者によるMATLABメモ| 双線形補間(バイリニア補間)のためのウェイト(重み)
このエントリに付けられたタグ|MATLABGeometry線形補間双線形補間バイリニア補間bilinear_weight重みウェイト画像

[基本] 連立方程式(線形)の解法

あまりに基本すぎて1回も書いてなかったので、一応書いておきます。
MATLABでの線形な連立方程式の解の求め方です。

これは至って簡単でMATLABにおける基本の演算子である「¥」(バックスラッシュ)を用いて実現できます。
簡単な例として下記の連立方程式を解いてみます。


この方程式を行列で表します。(ここでなんで?という人は線形代数の教科書とかを読んでください…)


このとき行列Aが逆行列を持つなら,未知の値x1、x2、x3は下記の行列の計算で求めることが出来ます。


実際にこの式をMATLABで書くと下記になり、解はx1=7、x2=-2、x3=-6であることがわかります。
x=[1,1,1];[0,2,-1];[1,0,1]]¥[-1;2;1]
x =
     7
    -2
    -6

なにかいい例を見つかりませんが、例えば下記参考記事内でベジエ曲線の制御点を求める際にもう少し複雑な線形連立方程式を解いています。
式の行列の作り方など多少参考になるかもしれません。

参考:MATLAB:マニュアル:MATLAB:関数:サイバネットシステム
参考:初心者によるMATLABメモ| 通過点列を通るC2連続な3次ベジエ曲線の生成
このエントリに付けられたタグ|MATLAB方程式連立方程式線形方程式解法

[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] 非線形方程式を解く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ヤコビ行列ヤコビアン行列

[Geometry] 通過点列を通るC2連続な3次ベジエ曲線の生成

指定した点列の間をC2連続に接続したベジェ曲線で補間し、そのベジェ曲線の制御点を返す関数を紹介します。
生成されるスプライン曲線は,C2連続ですので,結果的にG2連続(曲率連続)になっています。

簡単な使い方,やっていること、使用例の順に示します。(ダウンロードはこちらから)

まず、関数の使い方だけ書いておきます。
CPS = find_c2_bezier3(P) … 点列Pを厳密に通るC2連続な3次ベジエ・スプラインの制御点をセル配列CPSを返します。

以下に,まず簡単にやっていることの説明を書いておきます。
通過点p={P0...Pn}があるとき,点Piと点Pi+1の間を結ぶ3次ベジエ曲線の制御点を下記とします。
  (1)

3次ベジエ曲線をC2連続に接続するには,2次微分が連続となるので下記の条件(2)を満たせばいいことになります。
  (2)

式(2)は,3次ベジエ曲線の性質より,下記の式(3)となります。
  (3)

この式を未知数を含む項を左辺に,既知の値を右辺に寄せます。
  (4)

また、このままでは解と未知数の数が合いませんので,ここでは両端の2次微分が0になるという条件を追加しています。
  (5)

こうすることで解と未知数の数があい,線形な連立方程式として,未知の制御点ai,biを求めることができます。

以下が,この計算を行うMATLABの関数find_c2_bezier3()です。
上記の連立方程式を太文字の部分で解き,制御点を求めています。
function cps = find_c2_bezier3(pts) 

sz = size(pts);

dim = sz(1);
cnum = sz(2) - 1;
unnum = 2*cnum*dim;

%------ 連立方程式を立てる
eql = zeros(unnum);
eqr = zeros(unnum,1);

qq = reshape([1:unnum],dim,2*cnum);
qa = qq(:,1:2:end);
qb = qq(:,2:2:end);

% 始点の処理 C2 
for dg=1:dim 
    eql(dg, qa(dg,1)) = -2;
    eql(dg, qb(dg,1)) = 1;
end;
eqr(1:dim) = -pts(:,1);

for ii=1:cnum-1 
    jj = (2*ii-1)*dim;

    for dg=1:dim 
        % C1 condition
            eql(jj+dg, qb(dg,ii)) = 1;
            eql(jj+dg, qa(dg,ii+1)) = 1;

        % C2 condition
            eql(jj+dg+dim, qa(dg,ii)) = 1;
            eql(jj+dg+dim, qb(dg,ii)) = -2;
            eql(jj+dg+dim, qa(dg,ii+1)) = 2;
            eql(jj+dg+dim, qb(dg,ii+1)) = -1;
    end;

    eqr(jj+[1:2*dim]) = [2*pts(:,ii+1);zeros(dim,1)];
end;

% 終点の処理 C2
for dg=1:dim 
    eql((2*cnum-1)*dim+dg, qa(dg,end)) = 1;
    eql((2*cnum-1)*dim+dg, qb(dg,end)) = -2;
end;
eqr(end-dim+1:end) = -pts(:,end);

%------ 連立方程式を解き,制御点を求める
cps0 = eql¥eqr;
cps0 = reshape(cps0,dim,numel(cps0)/dim);

% 制御点をセグメントごとに整形
cps = cell(1,cnum);
for ii=1:cnum 
    cps{ii} = [pts(:,ii), cps0(:,2*(ii-1)+[1:2]), pts(:,ii+1)];
end;

以下,下記のような適当な点列をC2連続なベジエ曲線で補間してみます。
>> pts = [[0;0;0],[1;1;0.5],[2;1.5;2],[4;2;1.5],[5;2;-1],[5.5;0;0.5],[6.5;1.0;1.2]]
pts =
         0    1.0000    2.0000    4.0000    5.0000    5.5000    6.5000
         0    1.0000    1.5000    2.0000    2.0000         0    1.0000
         0    0.5000    2.0000    1.5000   -1.0000    0.5000    1.2000



この点列を補間し,制御点を取得します。
以下のように,セル配列の要素ひとつが1つの3次ベジエ曲線の制御点になります。
>> cps=find_c2_bezier3(pts)
cps = 
    [3x4 double]    [3x4 double]    [3x4 double]    [3x4 double]    [3x4 double]    [3x4 double]

>> cps{1}
ans =
         0    0.3626    0.7252    1.0000
         0    0.3765    0.7530    1.0000
         0    0.0495    0.0990    0.5000

実際に曲線を表示してみるとこんな感じです。


曲率をプロットしてみるとこのように曲率が連続になっていることがわかります。


以下、曲線表示のコード。
ここで使われているbezier()という関数は,この記事に載せている関数です。
nn = length(cps) + 1;

divn=50;
tl=0:1/divn:1;
step = length(tl);
tt = zeros(1,step*(nn-1));
coords = zeros(3,step*(nn-1));

thrpoint = pts;

for m=1:nn-1
    coords(:,step*(m-1)+1:step*m)=bezier(tl,cps{m});
end;

% 曲線の描画
col=[0,0,1];
col2=[1,0,0];

set(0, 'defaultAxesFontSize', 10, 'defaultTextFontSize', 10);

figure( 'InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [2 300 440 320]);
xlabel({'x'}, 'FontName','Times New Roman');
ylabel({'y'}, 'FontName','Times New Roman');
zlabel({'z'}, 'FontName','Times New Roman');
view([-36 30]);
box('on');
hold('all');

gr1 = plot3(coords(1,:),coords(2,:),coords(3,:), 'LineWidth',1.3, 'DisplayName', 'Result Curve', 'Color',col);
gr2 = plot3(pts(1,:), pts(2,:), pts(3,:),'DisplayName','Through Point', ...
    'MarkerFaceColor',col2, 'MarkerEdgeColor',col2, 'Marker','o', 'LineStyle','none', 'Color',col2);

axis equal;
set(gca, 'PlotBoxAspectRatioMode', 'manual', 'DataAspectRatio',[1 1 1]);


参考:初心者によるMATLABメモ| ベジェ曲線とバーンスタイン基底関数
参考:初心者によるMATLABメモ| 連立方程式(線形)の解法
ダウンロード:このエントリのmファイル(find_c2_bezier3.m)
このエントリに付けられたタグ|MATLABGeometry幾何ベジエ曲線ベジェ曲線Bezier補間線形方程式C2連続G2連続

[数値解析] ヤコビ行列(数値、Fast Code)

指定した関数のヤコビ行列を数値的に求める関数jacobi_numeric_fast()です。
先に書いておきますが,Fast-codeなので,複雑な問題においてこのヤコビアン行列を用いた場合には収束率が落ちる,解の精度が悪くなるなどの問題が起こる可能性があります。
ただ,意外とヤコビアン行列を手計算するのも大変なので,それなりに有効な手段にはなるかと思います。

以下,簡単な使い方,関数、使用例の順に示します。(ダウンロードはこちらから)

まず、関数の使い方です。
J = jacobi_numeric_fast(FUNC. X) … Xのときの関数FUNCの数値的に求めたヤコビアン行列をJとして取得します。
J = jacobi_numeric_fast(FUNC. X, H) … 微小幅Hを用いて,Xのときの関数FUNCの数値的に求めたヤコビアン行列をJとして取得します。
J = jacobi_numeric_fast(FUNC. X, H, A, B) … 関数FUNCに追加の引数A、Bを与えてヤコビアン行列を計算します。

※FUNCは関数ハンドル(@関数名)です。

コードは以下です。
基本的には,単に微小幅hを用いて,前後の点で関数を求めて差分を取っているだけです。
function J = jacobi_numeric_fast(func, x0, hh, varargin) 

% 可変の引数を関数に渡すようにするチェック
func = fcnchk(func,length(varargin));

nn = numel(x0);

if nargin<3 | isempty(hh) 
    hh=(1e-6)*max(x0(:))*ones(1,nn);
elseif numel(hh)==1 
    hh=hh*ones(1,nn);
elseif numel(hh)~=nn 
    error('変数の変化幅が不正です。');
end;

J = zeros(nn);
dx0 = eye(nn)/2;

for ii=1:nn 
    dxi = hh(ii)*dx0(:,ii);
    J(:,ii) = (func(x0 + dxi,varargin{:}) - func(x0 - dxi,varargin{:}))/hh(ii);
end;

以下,使用例。
こちらのエントリと同じ方程式を,jacobi_numeric_fast()で数値的に求めたヤコビアン行列を使って解いてみます。
fsolve()のための方程式の関数は,下記として定義しました。
ヤコビ行列Jは下記の太文字の部分のようにjacobi_numeric_fast()で求めました。
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(:);

function [F,J] = funcsjn(x) 
    F =funcs(x);
    if nargout>1 
        J = jacobi_numeric_fast(@funcs,x);
    end;

以下,fsolve()の実行結果です。これくらい簡単な問題だと,収束はあまり遅くならないようです。
>> 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,jacobi] = fsolve(@funcsjn, 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.41971e-019   2.68786e-005      2.34e-009             2.5
     5          6    4.93038e-032   2.53541e-010      4.13e-016             2.5
Optimization terminated: first-order optimality is less than options.TolFun.

こちらのエントリで求めているように解x=[0.7331, 1.8608, 1.3641]ですので,このように今回jacobi_numeric_fast()で数値的に求めたヤコビ行列を用いても同じ解を得られていることが確認できます。
>> x
x =

    0.7331
    1.8608
    1.3641

一応最終的なヤコビ行列は下記になります。
>> jacobi
jacobi =
    1.4662    3.7216         0
    1.8608    0.7331   -1.0000
    1.3641         0    0.7331

このような簡単な例では,うまく収束しましたが,実際に複雑な問題においてこのヤコビアン行列を用いた場合には収束率が落ちる,解の精度が悪くなるなどの問題が起こる可能性があります。
これを頭の隅に置いておいて貰えると良いと思います。

参考:初心者によるMATLABメモ| 非線形方程式を解く2 (fsolve)
ダウンロード:このエントリのmファイル(jacobi_numeric_fast.m)
このエントリに付けられたタグ|MATLAB数値解析ヤコビアン行列ヤコビ行列Fast-code

[その他] 初心者向けのおすすめの1冊:最新 使える!MATLAB

MATLABの初心者向けの本については,いろいろ出ていてホントに初心者が買う場合は迷うかと思います。
どれも一長一短ですし、周りの人に聞いても初心者にはどれがよいかは言うことが違うのでなんとも難しいですが,
読んだ感じ「初心者には」いいなぁと思った本を,何冊か紹介したていきたいと思います。気が向いたときになってしまいますが…

とりあえず、僕がいまだに時々参考にしているのは,この本ですね。
 Amazon.co.jp: 最新 使える!MATLAB: 本: 青山 貴伸,森口 肇,蔵本 一峰

目次
第1章 MATLABの基本的な操作 … はじめて使う人のために
第2章 グラフィックス … データの視覚化
第3章 スクリプトファイル … M‐ファイル
第4章 微分・積分
第5章 微分方程式
第6章 データ解析
第7章 Simulink
第8章 制御理論(古典制御)への適用
第9章 C MEX‐ファイル
第10章 COMサポートとDDEインターフェース―データリンクを便利に使う

内容的には上に書いた目次の通りですが,もう少し内容を触れてみるとこんな感じです。
第1章は、本当に初心者の人向けのMATLABを動かしてみよう的な内容。
第2章は、ざっくり言ってグラフの描き方。
第3章は,独自のプログラムファイルの書き方、独自の関数の作り方など。
第4章、第5章は,まぁそのままです。簡単な例を示してMATLABを用いた微分・積分のやり方。微分方程式の解き方。
第6章は、ちょっと目次のタイトルからはわかりづらいですが,データの読み込み,書き込みなど。
とんで第9章は、MATLABからのC言語系のプログラムの呼び出し方などをざっくりと簡単に紹介。
第10章は,VBからMATLABを操作,ActiveXの使用などをこちらもかなりざっくりと紹介。
第7,8章については,僕は基本的にSimulinkは使わないのでその部分は読んでいないので,7,8章に興味がある人には僕の感想は参考にならないかと思います。

こんな感じで内容は,全般に渡っています。(ただ深くはない)
また,レイアウトも結構よいです。
他のプログラム言語でプログラムをしたことがある人なら,これを読めばだいたい基本的なことは出来るようになると思います。
使いこなせるようにはならないと思いますが…
ただし,例えば「forループって何?」とかいう人にはさすがにこれだけでは,キツイと思います。

また、この本も他の初心者向けの本同様,
「MATLABは動くようになったけど結局やりたいことができない…」という状況に陥るような気がします…
ある程度使えるようになった人には,やはりこちらで紹介している本がかなりおすすめです。
実際いまだに頻繁に参考にしていますし。
このエントリに付けられたタグ|MATLAB初心者本

[数値解析] 非線形方程式を解く(ニュートン・ラプソン法)

[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以下なら解とみなす。) 注意!!
TolXxの変化に関する終了許容範囲 注意!!

次に出力として与えられる計算アルゴリズムが停止した理由を示す整数exitflagの説明を下記に書いておきます。
1関数が解 xに収束したことを示します。
2解xの変化が指定した許容値より小さくなったことを示します。
3関数の変化が指定した許容値より小さくなった(あまり変化しなくなった)ことを示します。
0関数計算の回数、あるいは繰り返し回数が最大回数に達していることを示します。
-1解が発散もしくは振動していることを示します。

また,最適化に関する情報を含む構造体outputの説明は下記です。
iterations実行した繰り返し回数
funcCount関数の評価回数
message終了メッセージ
answer最終的な解x
Func目的関数の最終的な大きさ
maxFunc目的関数の要素のうち絶対が最大のものの大きさ
TolFun関数値に関する終了許容範囲(maxFuncがTolFun以下なら解とみなす。)
TolXxの変化に関する終了許容範囲

以下に実際に例題を示します。
ここでは、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ファイル

[グラフ、プロット] 棒グラフの色を1本ずつ(だけ)変更する

棒グラフの色を1本ずつ(もしくは1本だけ)変更する方法。
サイバネットのテクニカルFAQにも書いてあるんだけど,テクニカルFAQに書いてある方法のようにbar()を繰り返して呼び出す方法よりこっちの方法が多分楽だと思います。

結論から先に書くと、単純な棒グラフの色を変えたいのであれば,別にbar()を何回も繰り返さないでもこれ↓で十分なはずです。この方が何本も書く場合にも面倒くさくないですし。
bar(diag(data),'stack');

ポイント的には,通常の棒グラフをわざわざ積み上げ棒グラフとして表現するところにあります。
積み上げ棒グラフbar(data,'stack')は下記の2つの性質を持っています。
 ・それぞれの項目の要素(x軸に振られているものを項目と言ってます)が異なる色で表現される。
 ・入力データは,行が項目,列が各要素としたマトリクス形式で表現される。

そこで,元のデータをdiag()を使って対角行列にしてやることで,各項目が異なる要素のみに0でないデータを持つ積み上げ棒グラフ用のデータを作成します。
あとは,そのデータで積み上げ棒グラフを描いてやれば,あたかもそれぞれの色が違う棒グラフになるわけです。
ちょっとわかりづらい説明ですいません。
これは実際に、例を見てもらった方が早いと思います。

下記に簡単に例を示します。
データとしてはこんな感じの適当なデータを使いました。
>> tcks={'1月';'2月';'3月';'4月';'5月';'6月'};
>> data = 100*rand(1,length(tcks))

data =
   79.2207   95.9492   65.5741    3.5712   84.9129   93.3993

まず、単純にbar()で棒グラフを書いてみます。
こうするとご覧のとおり全部の棒の色が同じな棒で棒グラフが描かれます。
figure('InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [200 100 320 160]);
hold all;
bar(data);
set(gca, 'Title',text( 'String', '縦の棒グラフ', 'FontSize', 10, 'FontName', 'MSゴシック'), 'XTickLabel',tcks ,'XTick',1:length(tcks) ,'XLim',[0,length(tcks)+1]);



これに対して,上に書いた方法で棒グラフを描いてみます。
すると、ご覧の通り無事にそれぞれの棒の色が異なる棒グラフが描けます。
ここではカラーマップをjet()にしていますが,別にそのままでも構いません。
figure('InvertHardcopy', 'off', 'Color', [1 1 1], 'Position', [200 100 320 160]);
hold all;
bar(diag(data),'stack');
colormap(jet);
set(gca, 'Title',text( 'String', '縦の棒グラフ', 'FontSize', 10, 'FontName', 'MSゴシック'), 'XTickLabel',tcks ,'XTick',1:length(tcks) ,'XLim',[0,length(tcks)+1]);



上を見て気づいたと思いますが,この棒グラフの色はカラーマップで制御されているので,カラーマップを明示的に作成して指定してやれば任意の色にすることができるわけです。
以下はこの記事のez_gradation()を用いて,黄色から緑へグラデーションするカラーマップを作成して,指定した例です。
>> colors = ez_gradation([[1,1,0];[0,0.66,0]], numel(data))
colors =
    1.0000    1.0000         0
    0.8000    0.9320         0
    0.6000    0.8640         0
    0.4000    0.7960         0
    0.2000    0.7280         0
         0    0.6600         0

>> colormap(colors);



また,1個だけ(ここでは3本目だけ)色を変えたいのであれば,こんな感じです。
colors = repmat([0,0,1],numel(data),1);
colors(3,:) = [0,1,0];
colormap(colors);




参考:MATLAB:マニュアル:MATLAB:関数:サイバネットシステム
参考:MATLAB:テクニカルFAQ bar関数で描画される棒グラフを、1本ずつ任意の色に変更するにはどうすればよいですか?:サイバネットシステム
参考:初心者によるMATLABメモ| カラーベクトルから手動でグラデーションを作る
このエントリに付けられたタグ|MATLABグラフプロット棒グラフカラーマップ色変更

[雑記,愚痴] 表示記事件数を変更しました

ちょっと意外と各記事が長くて1ページに10個の記事を表示するとあまりに1ページが長くなるので、 表示記事件数を5件に変更しました。
検索エンジンから来た人は,みたいページがみれないかもしれません。
その場合は,お手数ですが横のメニューのブログ内検索で検索してください。
このエントリに付けられたタグ|表示記事件数を変更

[基本] tic、tocで表示単位を指定して表示する関数

MATLABで処理の計算時間を測定するときには,ストップウォッチ関数tic()とtoc()を使っている人が多いと思います。
しかし,toc関数は下記のように秒単位での表示しか出来ません。
>> toc
経過時間は5353.697833秒です

軽い計算しかしない人は別に構わないのかもしれませんが,自分も含め数十分,数時間などの処理をする人の場合これでは困ります。
そこで,tocの表示形式を時間、分などで表示できるようにした関数tocf()を紹介します。
ダウンロードしてパスの通っているところに置いておくと,意外と便利だと思います。(ダウンロードは一番下のリンクから)
ダウンロードしてもらえればわかりますが,コードは別にたいしたことをしていないのでここでは紹介しません。

使い方は,「toc」と書いてあるところをそのまま「tocf」にするだけでO.K.です。
tic;
...
(メインの計算)
...
toc; → tocf;

するとこのような感じで秒だけでなく,時間,分などが表示されるようになります。
>> tocf
経過時間は1時間29分0.590977秒です

引数として整数値を与えると下記の表に対応した表示形式で表示されます。
整数値 表示形式
1経過時間は11.123456秒です
2経過時間は12分11.123456秒です
3経過時間は1時間12分11.123456秒です
4経過時間は1日1時間12分11.123456秒です
5Elapsed time is 12.123456 seconds.
6Elapsed time is 12 min 11.123456 sec.
7Elapsed time is 1 hr 12 min 11.123456 sec.
8Elapsed time is 1 day 1 hr 12 min 11.123456 sec.

実際に整数値を指定してやると,例えば分までの表示や英語表示ができるようになります。
>> tocf(2)
経過時間は89分4.148864秒です

>> tocf(6)
Elapsed time is 89 min 20.217385 sec.

また,引数が1のときはtocと同じになります。
>> tocf(1)
経過時間は5353.697833秒です

さらに下記のように出力を要求すると,経過日数、時間、分、秒の要素を持つ数値配列が得られます。
これを使って好きな形式で表示することも可能です。
>> tm=tocf
tm =
         0    2.0000   41.0000   48.8068
>> sprintf( '%d%s%d%s%02d%s%04.6f%s', tm(1), '日', tm(2), '時間', tm(3), '分', tm(4), '秒')
ans =
0日2時間41分48.806842秒


参考:MATLAB Function Reference - tic, toc
ダウンロード:このエントリのmファイル 左クリックでお願いします。
このエントリに付けられたタグ|MATLAB時間計算時間tictoctocf

[数値解析] 高次方程式(2次方程式)を解く

MATLABで1変数の高次方程式を解く方法ときには,roots()が使えます。
手始めに2次方程式で説明します。


MATLABでは,ここに書いてあるように,次数の高い項から行ベクトルとして係数を書くことで,多項式を表現できます。
つまり上記の2次方程式なら,下記です。
C=[a, b, c]

これを頭にいれて,下記の方程式を解いてみます。


この方程式は,次の形に変形できるので解は-1と-2のはずです。


実際に方程式をCとして表現し、roots()を利用して解いてると,下記のように,解は列ベクトルとして出力されます。
解も先に書いたように,-2と-1であることが確認できます。
>> C=[1 3 2];
>> roots(C)
ans =
    -2
    -1

また,重解の場合には下記のように2個同じものが解として得られます。
>> roots([1 2 1])
ans =
    -1
    -1

これが不都合な場合は,下記のようにunique()で除去してください。
>> unique(roots([1 2 1]))
ans =
    -1

また,roots()で解を求めると下記のように複素数の解も求まります。
>> A=roots([1 -2 4 4])
A =
   1.3425 + 2.0092i
   1.3425 - 2.0092i
  -0.6850  

この場合に実数解のみが必要な場合には下記のようにすればよいでしょう。
>> B=A(imag(A)==0)
B =
   -0.6850

また,コレと同様に,下記のような高次多項式で表される方程式についても同じように解くことが出来ます。


この多項式なら,下記のような形で表現します。


例えば,下記の5次方程式を解いてみます。


まず方程式は下記で表現されます。
>> C = [1, 8, 18, 4, -19, -12];

この場合もroots()を用いて下記のように解くことがます。
>> roots(C)
ans =
   -4.0000
   -3.0000
    1.0000
   -1.0000
   -1.0000


参考:MATLAB Function Reference - roots
参考:MATLAB:マニュアル:R14:MATLAB:数学:サイバネットシステム-多項式の表現
参考:MATLAB Function Reference - imag
このエントリに付けられたタグ|MATLAB数値解析方程式高次方程式二次方程式一元二次方程式roots

[基本] 行列(ベクトル)の大きさを調べる

プログラムを組んでいて,行列の大きさが知りたい場面は多いと思います。
そんなときに使う関数のまとめです。かなり初歩的なことなので、あまり書くこともありませんが。
MATLABでは,下記の4種類を使うことによって行列の形状と大きさを取得することが出来ます。一応一番右の列の別の表現も覚えておくと便利だと思います。
個人的にはlength()は普通のプログラミング言語のそれと異なる挙動なので注意が必要かな?と思います。
関数名関数の説明一般的に等価な処理(別の表現)
size()各次元ごとの要素数を配列として返します。 
length()行の各次元の要素数のうち最大のものを返す( 0 のサイズの行列または配列の長さは 0)。
列ベクトル,行ベクトルに対してはnumel()と同じ。
max(size(M))と等価。
ndims()行列の次元を取得。(最低2次元)length(size(M))と等価。
numel()行列の全要素数を取得。列ベクトル,行ベクトルに対してはlength()と同じ。length(M(:))と等価。またprod(size(M))とも等価。

以下,下記の3×4の行列Aを使って簡単な例を示します。
>> A=floor(10*rand(3,4))
A =
     5     5     4     2
     6     4     8     2
     4     1     8     5

まず、size()を使ってこの行列の各次元の要素数を調べます。
>> size(A)
ans =
     3     4

続いて,length()を使って,各次元の要素数のうち最大のものを取得します。
>> length(A)
ans =
     4

以下は,length()の別の表現。
>> max(size(A))
ans =
     4

続いて,numel()を使って,行列の全要素数を調べます。
>> numel(A)
ans =
    12

以下の2つは,numel()の別の表現です。
>> prod(size(A))
ans =
    12

>> length(A(:))
ans =
    12

最後に、ndims()を使って行列が何次元かを調べます。この場合は3×4なので2次元行列であることがわかります。
>> ndims(A)
ans =
     2

ndims()は,下記の表現と同じです。
>> numel(size(A))
ans =
     2

>> length(size(A))
ans =
     2

また,例えば一列の長さ(2次元行列の縦の長さ)を調べたいなら、下記です。
>> size(A,1)
ans =
     3

>> length(A(:,1))
ans =
     3

同様に,一行の長さ(2次元行列の横の長さ)を調べたいなら、下記です。
>> size(A,2)
ans =
     4

>> length(A(1,:))
ans =
     4

このように,size(A,n)とすることで調べたい次元を指定して要素数を調べることが出来ます。
これは多次元行列の場合も同じです。

参考:MATLAB Function Reference - size
参考:MATLAB Function Reference - ndims
参考:MATLAB Function Reference - numel
参考:MATLAB Function Reference - length
このエントリに付けられたタグ|MATLAB行列の大きさ大きさ長さサイズsize()numel()length()ndims()

[雑記,愚痴] MATLABのメニューバーを見ていて思った

そういえばMATLABってGUIで意外と出来る操作も意外とあるなぁ。
Figure弄るところぐらいしか使っていないや、自分。
まぁGUI使っちゃうと結局バッチ処理できなくて面倒だから,今後も本格的に使うことはないんだろうけど。
ちょっとしたことやるだけなら、覚えた方がいいのかも…
このエントリに付けられたタグ|雑記

[基本] 起動時にデフォルトカラーなどを設定する (起動時にスクリプトを実行)

MATLABでは,MATLABのホームディレクトリにstartup.mファイルが存在すれば、起動する際にstartup.mが実行されます。
これを利用することで,毎回パスを追加しなおしたりする手間が省けます。
他にもデフォルトのフォントやサイズを変更したり,自分が最近メインで使っているワークディレクトリに移動したりすることができます。

一応サイバネットのサイトを見ると,startup.mを置く場所は下記だそうです。(実際にはパスが通ってればこれに限りませんが)
それとWindowsでMy Document以下にmatlab
[Windows版の場合]
ショートカットで指定しているMATLABの作業フォルダ。MATLABがデフォルトで自動的に作成する作業フォルダは、次の通りです。
* R2007a:各ユーザのMy Documentフォルダの下のMATLABフォルダ
* R2006b以前:MATLABインストールフォルダの下のworkフォルダ

[非Windows版の場合]
自動的にMATLABパスに追加される$HOME/matlabディレクトリ下、もしくは環境変数MATLABPATHで定義したディレクトリ。
参考:MATLAB:技術サポート:テクニカルFAQユーザ毎にMATLABの環境をカスタマイズ出来ますか?:サイバネットシステム

単純にmファイルを実行するだけなので,ある意味なんでも設定できますが、メニューバーの「ファイル(File)→設定(Preferences)」から起動できる設定ダイアログで設定できるものはそちらで設定して,startup.mにはそれ以外の設定を書いておくのがいいかなというのが私見です。
で、例えば自分用のライブラリフォルダなどは,日々フォルダが作成されたり,階層が掘られたりするのでstartup.mからパスを追加してやった方が楽だと思います。
また、Figureについての設定は結局個別に弄るでしょうから気休め程度です。

以下,一応startup.mの例を載せておいてみます。
書きすぎた気もしますし、これをそのまま使うとややFigureは変になると思います。
僕自身は最初の3つと多少のFigure設定(Figureウインドウの背景を透明になど)しか実際にはやっていません。
これ以外にももろもろの設定を行うことができます。
%%%% 作業フォルダに移動
cd('C:¥sach1o¥MATLAB');

%%%% 自分の作ったライブラリを全部追加(指定したディレクトリ以下を全部パスに追加)
addpath(genpath('C:¥sach1o¥MATLAB¥lib'));

%% diaryコマンドが作るファイルに拡張子を追加
set(0, 'DiaryFile', 'diary.txt');


%%%% Figureウインドウ関係のデフォルト設定 -- 以下は参考程度に

scrsz = get(0, 'ScreenSize');    % スクリーンサイズの取得

% Figureウインドウの位置とサイズのデフォルト
set(0, 'defaultFigurePosition', [200 scrsz(4)-620 480 360]);

% Figureの紙の大きさ
set(0, 'defaultFigurePaperType', 'A4');

% Figureの寸法の単位
set(0, 'defaultFigurePaperUnits', 'inches');

% 黒背景のものprintするときに色を変更するか?
set(0, 'defaultFigureInvertHardcopy', 'on');

% 軸のボックスをonに
set(0, 'defaultAxesBox', 'on');

%% 色関係
% テキストカラー
set(0, 'defaultTextColor', [1 0 0]);

% 各軸の色 (軸の目盛りのテキストも変更)
set(0, 'defaultAxesXColor', [0 0 0]);
set(0, 'defaultAxesYColor', [0 0 0]);
set(0, 'defaultAxesZColor', [0 0 0]);

% パッチのエッジの色
set(0, 'defaultPatchEdgeColor', [0 0 0]);

% Surfaceのエッジの色
set(0, 'defaultSurfaceEdgeColor', [0 0 0]);

% 注釈オブジェクトの線の色 plot の色ではなくannotationの色 Figureウインドウで書ける線の色
set(0, 'defaultLineColor', [1 0 0]);

% Figureウインドウの色
%    この例では透明にしていますが、[1 1 1]なら白などカラーベクトルで指定できます。
set(0, 'defaultFigureColor', 'none');

% 軸の色(プロットエリアの色)
set(0, 'defaultAxesColor', [1 1 1]);

% plotなどで自動的に付けられる色の指定
corder = [[1 0 0];[0 0.5000 0];[0 0 1];[1 0.7500 0.7500];[0.7500 0 0.7500];[0.7500 0.7500 0];[0.2500 0.2500 0.2500]];
set(0, 'defaultAxesColorOrder', corder);

% デフォルトのカラーマップ
cmap = spring(128);
set(0, 'defaultFigureColormap', cmap);

%% フォント関係
% GUIのフォント
set(0, 'defaultUicontrolFontName', 'MS UI Gothic');
% 軸のフォント
set(0, 'defaultAxesFontName', 'Arial');
% タイトル、注釈などのフォント
set(0, 'defaultTextFontName', 'Times');

% GUIのフォントサイズ
set(0, 'defaultUicontrolFontSize', 9);

% 軸のフォントサイズ
set(0, 'defaultAxesFontSize', 10);

% タイトル、注釈などのフォントサイズ
set(0, 'defaultTextFontSize', 12);

%% そのほか
% 軸の線の太さ
set(0, 'DefaultAxesLineWidth', 3);

% 注釈オブジェクトの線の太さ
set(0, 'DefaultLineLineWidth', 3);

また、起動時にオプションを指定してやる方法(これとかこれ)もありますが、こちらはMATLABを他の言語から使用したい人以外はあまり使わないと思います。

参考:MATLAB:マニュアル:MATLAB:関数:サイバネットシステム
参考:MATLAB:技術サポート:テクニカルFAQユーザ毎にMATLABの環境をカスタマイズ出来ますか?:サイバネットシステム
参考:MATLAB:マニュアル:MATLAB:Desktop Tools and Development Environment:サイバネットシステム
参考:MATLAB:マニュアル:MATLAB:Desktop Tools and Development Environment:サイバネットシステム
このエントリに付けられたタグ|MATLABstartup()環境設定Figure

[小技] diary()を引数なしで使ったときのファイル名に拡張子をつける

先ほどstartup.mのエントリを書いていて思ったんですが,一部抜き出しでこちらに書いておきます。
diary()を使うと,コマンドウインドウの出力のログを作成してくれますが,この関数は引数なしで使ったときの作られるファイル名がなぜか拡張子なしの「diary」なんですよね。
ログはテキストファイルで出力されるので,拡張子が付いて欲しい…

その場合には,下記のコマンドをstartup.mに追加しておくと,拡張子.txtが付くようになります。
%% diaryコマンドが作る標準のファイル名を変更
set(0, 'DiaryFile', 'diary.txt');


参考:MATLAB Function Reference - diary
参考:初心者によるMATLABメモ| 起動時にデフォルトカラーなどを設定する (起動時にスクリプトを実行)
このエントリに付けられたタグ|MATLABdiary()デフォルトstartup()

[関数] 引数の数が可変な関数

引数が決まった数でない関数を作りたい場合に使うのは,vararginです。
変数vararginは,セル配列であり,これを利用することで関数に対して、任意の数の引数を与えることを可能にします。

以下、簡単な例です。
まず、簡単に2つの引数の和を計算する関数vararg1()を考えます。
function val = vararg1(in1, in2)
    val = in1 + in2;

>> vararg1(1,2)
ans =
     3

同様に3つの引数のの和を計算する関数vararg2()はこんな感じになります。
function val = vararg2(in1, in2, in3)
    val = in1 + in2 + in3;

>> vararg2(1,2,3)
ans =
     6

このように3つくらいの引数なら別段問題はないと思うのですが,例えば引数が1個から100個まで引数を持つ関数を作る場合にはどうするのか?
まぁさすがに100個は極端ですが,そのように引数の数が変わる可能性がある関数を作成するときに使うのが、vararginです。
例えば,すべての引数の総和を求める関数はvararginを利用して、下記のように書けます。
function val = vararg3(varargin)
    in = [varargin{:}];
    val = sum(in);

実際にこれを実行するとこんな感じになり、任意の数の引数を持つ関数になっていることが分かると思います。。
一応、この例では2行目のセミコロンを外してinの値を表示しています。
>> vararg3(1,2,3,4,5,6)
in =
     1     2     3     4     5     6

ans =
    21

また、同様に例えば1個目の引数だけは名前を最初から付けておきたい場合は下記のような関数の宣言にすることで実現できます。
この場合は一つ目の引数だけ、変数in1に格納されます。
function val = vararg4(in1, varargin)
    in = [varargin{:}];
    val = in1 + sum(in);

先程の例と同様に,2行目のセミコロンを外してinの値を表示してみると、下記のように第一引数はvararginに格納されていないことが確認できます。
>> vararg4(1,2,3,4,5,6)
in =
     2     3     4     5     6

ans =
    21

また、vararginは単にセル配列の型の変数なので、下記のように直にアクセスして値を変えたりすることもできます。
function val = vararg5(varargin)
    varargin{1} = 11;
    in = [varargin{:}];
    val = sum(in);

最後にこの表現は,下記のような無名関数の場合にも有効です。
vararg6 = @(varargin) sum([varargin{:}]);


参考:MATLAB Function Reference - varargin, varargout
参考:初心者によるMATLABメモ| 出力の数が可変な関数
このエントリに付けられたタグ|MATLAB関数引数可変引数varargin

[関数] 出力の数が可変な関数

引数が決まった数でない関数を作りたい場合に使うのはvarargin,でしたが、同様に出力が決まった数でない場合はvarargoutを使います。
変数varargoutは,vararginと同様にセル配列であり,これを利用することで任意の数の出力を関数が出すことをことを可能にします。

以下、簡単な例です。
まず、簡単に2つの出力(ここでは引数の1乗と2乗)を計算する関数func1()を考えます。
function [val1, val2] = func1(in)
    val1 = in;
    val2 = in.^2;

>> [v1, v2] = func1(2);
>> [v1, v2]
ans =
     2     4

同様に3乗も出力する関数func2()はこんな感じになります。
function [val1, val2, val3] = func2(in)
    val1 = in;
    val2 = in.^2;
    val3 = in.^3;

>> [v1, v2] = func2(2);
>> [v1, v2]
ans =
     2     4

>> [v1, v2, v3] = func2(2);
>> [v1, v2, v3]
ans =
     2     4     8

このように3つくらいの出力なら上記のように書き加えていけば問題ありません。
このように引数の数が変わる可能性がある関数を作成するときに使うのが、varargoutです。
例えば,出力の数にあわせて引数のn乗を出力する関数はvarargoutを利用して,下記のように書けます。
function varargout = func3(in)
    nn = max(nargout,1);
    varargout = num2cell(in.^[1:nn]);

実際にこれを実行するとこんな感じになり、任意の数の出力を持つ関数になっていることが分かると思います。
>> func3(2)
ans =
     2

>> [v1, v2] = func3(2);
>> [v1, v2]
ans =
     2     4

>> [v1, v2, v3, v4, v5] = func3(2);
>> [v1, v2, v3, v4, v5]
ans =
     2     4     8    16    32

実際にはfunc3のように書くより、下記のようにセル配列に直にアクセスする感じで書くことになると思います。
function varargout = func4(in)
    nn = max(nargout,1);
    for ii=1:nn 
        varargout{ii} = in.^ii;
    end;

>> [v1, v2, v3, v4, v5, v6] = func4(2);
>> [v1, v2, v3, v4, v5, v6]
ans =
     2     4     8    16    32    64

参考:MATLAB Function Reference - varargin, varargout
参考:初心者によるMATLABメモ| 引数の数が可変な関数
このエントリに付けられたタグ|MATLAB関数出力可変出力varargout