Links:
Matlabでは直接アニメーションgifは作れないようなので,一旦tif形式で1枚1枚のファイルに書き出して, image magic付属の(gradsにもついている)convertでファイルの結合・gifへの変換・アニメーション化,を行う. サンプルプログラムを示そう.
tsz=100 dim=ones(20,10); xs=[1:1:20]'*ones(1,10); ys=ones(20,1)*[1:10]; k=2*pi/10; l=2*pi/5; omega=2*pi/10; for t=1:tsz dim=sin(k * xs + l * ys - omega * t); tifffile=sprintf('d:\\tmp\\%02d.tif',t) contour(dim',[-0.8 -0.6 -0.4 -0.2 0.2 0.4 0.6 0.8]) print(tifffile,'-dtiff','-r75') hold off; end cd d:\tmp pwd pause !convert -loop 0 -delay 20 *.tif anim.gif % loopを指定しないとpower pointでanimationしない.
海洋の深度や大気の気圧は,鉛直上向きに減少ししかもマイナス記号をつけないのが一般的だ.これをMatlabの3D Graphで実現するには,
1) 妥協してマイナス記号をつける
2) 図の下向きが,物理空間の鉛直上向きで一旦図を描き,Figureをクリックしてz軸のプロパティの「降順」をクリックする.
の2通りがある.(1)の方が安直なので自分で見るには適しているが,最終的な図を作るには(2)がよい.
以下を実行して見て取れるように,マニュアルには載っていないけれど(と思う),線の太さなどもplotコマンドで指定できる. Matlabの線はデフォルトでは細すぎて,よく見えないことが多いので,適宜太くするとよいだろう.
plot(t,sin(2*t),'-mo',... 'LineWidth',2,... 'MarkerEdgeColor','k',... 'MarkerFaceColor',[.49 1 .63],... 'MarkerSize',12)
Matlabでは行列の計算は定義されているが,テンソルになるとダメである.例えば,3次元配列(2×2×2)と1次元配列の掛け算はやってくれない(エラーになる).もっと簡単な,
x=ones(2,2,2); y=ones(2,1); z=x(1,1,:)*y;つまり(1×1×2)と(2×1)をやらせようとしても
??? エラー: ==> * 行列入力はサポートしていません.となってしまう.このような場合は,reshapeを使って,行列の演算にしてやることで計算ができる.つまり,
z=reshape(x(1,1,:),1,2)*y;とすれば良い.
sscanfで空白を含む文字を読み出すと、%cであれば空白が保存されるが、%sでは除去される。 通常空白を保存したいだろうから、%c(%sではなく)を使う方が良い。書き出しはどっちでも良い。次の例で動作の違いが分かるだろう。
cstr='a b c'; sscanf(cstr,'%s') sscanf(cstr,'%c') sprintf('%s',cstr) sprintf('%c',cstr) ans = abc ans = a b c ans = a b c ans = a b c
Matlabの欠損値 NaN を含むデータを、他のソフトウェアで利用するには(例えばGradsでの描画)、NaN を適当なユーザーが指定した定数に変換してやらなくてはならないことがしばしばある。この変換は、次のようにしてmatlab の isnan関数を使うことで可能である。
x=[1 NaN 3]; x(isnan(x))=0 x = 1 0 3
ちなみにmatlabではNaNとの==,~=などの評価は機能しない。たとえば、
x=[1 NaN]; x==NaN x~=NaN ans = 0 0 ans = 1 1
のように、NaN との==による比較は常に偽、~=での比較は常に真が返される。したがって、 NaN以外の値であれば、x(x==3)=0 などとして変換ができるが、NaNについてはそうはいかずisnanを使わざるを得ないことを記憶しておくべきである。なお、無限大になる可能性のある演算を含んでいる場合は、NaNとともにinfを除外するために、isfinite関数を使うこともできる。ただし、無限大になる可能性のある演算は本来行うべきではないので、isnanをもちいる方がゼロ割の間違いを発見しやすい。
単純な転置'では、複素数は転置に加えて共役が取られる。転置だけしたい場合は、.'を用いること。例えば
x=[1+i 2+2*i; 3+3*i 4+4*i]; x' x.'
を実行すると
ans = 1.0000 - 1.0000i 3.0000 - 3.0000i 2.0000 - 2.0000i 4.0000 - 4.0000i ans = 1.0000 + 1.0000i 3.0000 + 3.0000i 2.0000 + 2.0000i 4.0000 + 4.0000iとなる。
xが列ベクトル、行ベクトルのどちらでも、以下の操作で必ず列ベクトルになる。
x = x(:);
ただし、(xsz,yxz)の行列も、(xsz*ysz,1)の列ベクトルに変換されるので、行列への 利用は慎重に。
[path,name,ext] = fileparts(outctlname) outctl.grdname=sprintf('^%s.grd',name)
関数の引渡し変数にundefを入れておいて,それが定義されていないことを次のようにしてチェックし,適宜な値を入れる.
if exist('undef')~=1; undef=-1.e+10; end
重要なのはメッセージを出してそこで処理を止めること.
たとえば、if (strcmpi(inctl.tstp(2:3),'MO') & strcmpi(inctl.tstp(2:3),'YR') '*testfun* unit of tstp is not MO nor YR',var,pause end
引数が整数の場合 | |
d | 引数を10進値として出力する |
u | 引数を符号無し10進値として出力する |
o | 引数を符号無し8進値として出力する |
x | 引数を符号無し16進値として出力する |
引数が浮動小数点の場合 | |
f | 引数をddd.ddd(dは10進数の1桁)の形式で出力する |
e | 引数をd.ddde+ddの書式で出力する |
A(1:nsiz,1)*ones(1,msiz)あるいは
ones(nsiz,1)*A(1,1:msiz)というようように、onesを用いて行列の掛け算をする方法と、
A(1:nsiz,ones(1,msiz))あるいは
A(ones(nsiz,1),1:msiz)のように行列の中にonesを書く方法とがある。一般に後者の 法が、速度が速いようである。ちなみに
A(ones(n1siz,n2siz),1:msiz) %n1siz,n2sizを別に指定しても A(ones(n1siz*n2siz,1),1:msiz) %n1siz*n2sizのみが意味あるなら、別々に指定するべきではないの二つで得られる結果は変わらない。最初の例は、プログラムの可読性を損なうので、行うべきではない。
多くの関数は第一要素について演算を行う。
第一要素について演算を行うことが、
size(sum(randn(3,2))) size(prod(randn(3,2))) size(mean(randn(3,2))) size(median(randn(3,2)))
の答えがすべて、[1 2] となることで確認できる。
x=[3 7 5; 0 4 2]; sort(x)
とすると答えが、
0 4 2 3 7 5
となることで確認できる。
第一要素について演算を行う。 fft のhelp には「FFTは、行列に対しては列単位で操作を行います。」となっている。 これは、dim(xsz,ysz)という配列があれば、xに関して(yを固定して)演算を行うということである。 これはtaper関数(hamming等)が[nsiz,1]の大きさの結果を返すことと整合的である。
% 列について変数 x=randn(1,8); x(2,:)=x; f=fft(x); max(abs(f(1,:)-f(2,:))) %ゼロならf(1,:)とf(2,:)は同一 % 行について変数 x=randn(8,1); x(:,2)=x; f=fft(x); max(abs(f(:,1)-f(:,2))) %ゼロならf(:,1)とf(:,2)は同一
関数中で
if (nargout==0) plot(x1,y1) end
のようにすると、出力引数が無い場合に画面に図を出す関数ができる。
/opt/matlab/toolbox/stats
岡田です。 まず、Matlab で 3D の図が書かれている tiff ファイル群を作るスクリプト (2つ)をお送りします。 ・図を描くのに必要なデータ データを収めていたマシンが修理から帰ってきたばかりなので、データを 引き出せません。今日の夜にはデータを出せると思います。 ・tiff ファイルから GIF Animation を作るプログラム。 は必要でしょうか。お知らせください。 ---------------------------------------------------------------------- スクリプト 1 %cd d:\kanamoto\matlab\bd200 clear for I=121:174 J=0; II=28 JJ=35 isz=20; i1=1+2*(II-1); i2=i1+isz-1; jsz=20; j1=1+2*(JJ-1); j2=j1+jsz-1; ksz=20; k1=1 ; k2=k1+ksz-1; imsh=2;jmsh=2;kmsh=2; cmax= 2.5E-04; casename='bd20b' % filename=sprintf('%s%03d','ut.',I); filename=sprintf('%s%s%03d',casename,'ut.',I); fid=fopen(filename,'r','ieee-be');rdsize=[128*128*20]; ur=fread(fid,rdsize,'float32'); ur=reshape(ur,128,128,20); fclose(fid); % filename=sprintf('%s%03d','vt.',I); filename=sprintf('%s%s%03d',casename,'vt.',I); fid=fopen(filename,'r','ieee-be'); vr=fread(fid,rdsize,'float32'); vr=reshape(vr,128,128,20); fclose(fid); % filename=sprintf('%s%03d','wt.',I); filename=sprintf('%s%s%03d',casename,'wt.',I); fid=fopen(filename,'r','ieee-be'); wr=fread(fid,rdsize,'float32'); wr=reshape(wr,128,128,20); fclose(fid); % filename=sprintf('%s%03d','o.',I) filename=sprintf('%s%s%03d',casename,'o.',I) fid=fopen(filename,'r','ieee-be'); omegar=fread(fid,rdsize,'float32'); omegar=reshape(omegar,128,128,20); fclose(fid); omegzr=zeros(size(ur)); omegzr=omegaz(ur,vr); for K=k1:k2 u(:,:,K)=ur(i1:i2,j1:j2,K)'; v(:,:,K)=vr(i1:i2,j1:j2,K)'; w(:,:,K)=wr(i1:i2,j1:j2,K)'; omega(:,:,K)=omegar(i1:i2,j1:j2,K)'; omegz(:,:,K)=omegzr(i1:i2,j1:j2,K)'; end clear *r omga=zeros(size(omega)); omgb=zeros(size(omega)); omga(omegz>0)=omega(omegz>0); omga(omegz<0)=-1; omgb(omegz<0)=omega(omegz<0); omgb(omegz>0)=-1; cmin= cmax; imsh1=i1*0.05;imsh=imsh*0.05;imsh2=i2*0.05; jmsh1=j1*0.05;jmsh=jmsh*0.05;jmsh2=j2*0.05; kmsh1=k1*0.05;kmsh=kmsh*0.05;kmsh2=k2*0.05; x=((i1:i2)-0.5)*0.05;y=((j1:j2)-0.5)*0.05;z=((k1:k2)-0.5)*0.05; clf p = patch(isosurface(x, y, z, omgb, cmin)); set(p, 'FaceColor', [0,0,1], 'EdgeColor', 'none'); isonormals(x,y,z,omgb, p) p = patch(isosurface(x, y, z, omga, cmax)); set(p, 'FaceColor', [1,0,0], 'EdgeColor', 'none'); isonormals(x,y,z,omga, p) daspect([1 1 1]); [cx cy cz] = meshgrid(imsh1:imsh:imsh2,jmsh1:jmsh:jmsh2,kmsh1:kmsh:kmsh2); h2=coneplot(x,y,z,u,v,w,cx,cy,cz,2); set(h2, 'FaceColor', 'green', 'EdgeColor', 'none'); % xlim([min(x) max(x)]); ylim([min(y) max(y)]); zlim([min(z) max(z)]); xlim([(i1-1.00000001)*0.05 i2*0.05]); ylim([(j1-1.00000001)*0.05 j2*0.05]); zlim([0.0 1.0]); grid on; for ID=1:1 % dtheta=15*(ID-1) % dphi=15*sin((ID-1)*pi/12) %% dtheta=37.5 %% dphi=60 view(3) camva(11); lighting gouraud; lit=[1/ID^2,1/ID^2,1/ID^2]; light('Position',campos,'Color',lit); % camorbit(dtheta,dphi); % stringx=sprintf('%s%d%s%s%+10.1E%s','\bf ',I,' hour ','渦度=',cmax,'s^-^1'); stringx=sprintf('%s%d%s','\bf ',I,' hours'); title(stringx,'Fontsize',12) xlabel(' x (km)'); ylabel(' y (km)'); zlabel(' z (km)'); % axis tight; box on; % camva(11); J=J+1; tifffile=sprintf('%s%03d%03d%s','bd20butca',I,J,'.tif') print(tifffile,'-dtiff','-r75') end % cla % close all % clear p % clear o* % clear end ---------------------------------------------------------------------- スクリプト 2 以下のスクリプトは omegaz.m という名前にしておく。 function omegz=omegaz(u,v); omegz=zeros(size(u)); omegz(2:end-1,2:end-1,:)=v(3:end,2:end-1,:)-v(1:end-2,2:end-1,:)-(u(2:end-1,3:end,:)-u(2:end-1,1:end-2,:)); omegz( 1 ,2:end-1,:)=v(2,2:end-1,:)-v( end ,2:end-1,:)-u( 1 ,3:end,:)+u( 1 ,1:end-2,:); omegz(end,2:end-1,:)=v(1,2:end-1,:)-v(end-1,2:end-1,:)-u(end,3:end,:)+u(end,1:end-2,:); omegz(2:end-1, 1 ,:)=v(3:end, 1 ,:)-v(1:end-2, 1 ,:)-u(2:end-1,2,:)+u(2:end-1, end ,:); omegz(2:end-1,end,:)=v(3:end,end,:)-v(1:end-2,end,:)-u(2:end-1,1,:)+u(2:end-1,end-1,:); omegz=omegz/100.0;
peak1=peaks; peak2=peak1./3.; subplot(2,1,1); cint=1; lvmn=floor(min(min(peak1))/cint)*cint ; lvmx=ceil(max(max(peak1))/cint)*cint; clevs=lvmn:cint:lvmx; contourf(peak1,clevs); colorbar; shading flat; subplot(2,1,2); contourf(peak2,clevs); colorbar; shading flat;Thank you very much in advance. Shoshiro Minobe minobe@neptune.sci.hokudai.ac.jp Graduate School of Science, Hokkaido University, Japan ---------------------回答 6 Mar 1998--------------------------------- Hello, I haven't found a perfect solution to your problem, but here is a possible workaround:
peak1=peaks; peak2=peak1./3.; % make a copy of the last column lastrow=peak2(:,end); %append four copies of the last column to the data peak2(:,end+1:end+4)=lastrow(:,ones(1,4)); peak2(end,end)=max(max(peak1)); subplot(2,1,1); cint=1; lvmn=floor(min(min(peak1))/cint)*cint ; lvmx=ceil(max(max(peak1))/cint)*cint; clevs=lvmn:cint:lvmx; contourf(peak1,clevs); colorbar; shading flat; subplot(2,1,2); [cout hand]=contourf(peak2,clevs); % get the XLIM property of the subplot xlims=get(gca,'xlim'); %hide the last four columns by changing the xlim property set(gca,'xlim',[xlims(1) xlims(2)-4]); colorbar; shading flat;
If you have any questions, please send email (日本語OK) to me or for general MATLAB support to techmat@cybernet.co.jp Regards, Rob Henson Cybernet Systems MATLAB Distributor in JAPAN rob@cybernet.co.jp