ここに説明があるので,ここではHowTo的に描いてみる.
インストール
まずはMATLABあるいはOctaveをインストールする.起動すると, こんな感じだ:
ファイル形式
読み込めるファイルは, 次のようなものである:
X座標 Y座標 関数値1 関数値2 ....
で, 複数行書く時には, Y座標を変更していき, 次にX座標を変更させた方が手数が少なくて済む.
ソフトウェア操作
起動してフォルダーを選択し, ファイルを読み込む. 例えば my.data が
0 0 1.0 0 1 2.0 0 2 5.0 1 0 1.0 1 1 2.0 1 2 5.0
まあこの場合は\(u(x,y)=1+y^2\)ということである.この場合
>> A=load("my.data") >> x=A(:,1) #これで1列目がxベクトル >> y=A(:,2) #これで2列目がxベクトル >> u=A(:,3) #これで3列目がuベクトル
いやいや, それらは2次元データなので reshape するわけね.
>> x=reshape(x,3,2) #y方向に3格子, x方向に2格子 x =0 1 0 1 0 1 >> y=reshape(y,3,2) y =0 0 1 1 2 2 >> u=reshape(u,3,2) u =1 1 2 2 5 5
なるほど. では描いてみる:
>> contour(x,y,u,10)
デフォルトでは, かな〜りダサい図を描く:
ツッコミどころ満載である:
- 文字がダサい.まるでコンピュータが描きそうな字だ.Timesくらい使えよ
- set(gca, "FontName", "Times New Roman");
- 文字が小さい.講義室の後ろで「みえませーん」叛乱が起こる
- set(gca,'fontsize',20) で訂正可能
- 線が細い.
- 線のラベルがないので関数値が不明
- contour(x,y,u,10,'ShowText','on')で一応つくが, 使い物にはならない
- XとYのスケールが異なる
- pbaspect([1 1 1])で訂正可能
いろいろいじるのが面倒である.その場合には, 保存して別のソフトで描くことも可能である.
いつも使っているのはこれ:
function [val lev] = savecontour(filename,cnt) fid=fopen(filename,'w'); npar = 0; ptr = 0; N = size(cnt,2); while(ptr < N) npar++; ptr++; lev{npar} = cnt(1,ptr); n = cnt(2,ptr); val{npar} = cnt(:,((ptr+1):(ptr+n)))'; for i=1:n fprintf(fid,"%f %f\n",cnt(1,ptr+i),cnt(2,ptr+i)); end fprintf(fid,"\n"); ptr+=n; endwhile; fclose(fid); endfunction;
これで
[ cnt lev ]=contourc(x,y,u); savecontour("mycont.dat",cnt)
こうすると mycont.dat に X Y座標が保存されるので, あとは好きなように描けば良い. たとえば
graph -Tps -C -x 0 1 -y 0 2 -nx -ny --line-width 0.0025 --frame-lin\
e-width 0.002 mycont.dat |ps2pdf -dEPSCrop - new.pdf
これでPDFファイルをゲットだ:
解を保存
収束後の解を保存しよう. 例えばデータを保存するC++コードは
#include <boost/filesystem.hpp> ... std::ofstream ofs; ... boost::filesystem::path folder=getenv("HOME"); boost::filesystem::current_path(folder/"tes_test"); ofs.open(filename); ... void write() { for(int i=0;i<U.size(0);i++) for(int j=0;j<U.size(1);j++) ofs << static_cast<double>(i)/(U.size(0)-1) << " " << static_cast<double>(j)/(U.size(1)-1) << " " << U(i,j) <<std::endl; }
こんな感じである. サンプルの図は100x100格子で計算したが, 収束が遅く, 前回とのずれが\(10^{-10}\)に達するまでに約30000イタレーションを要した.
こいつを作図するOctaveコードを, jacobi.m ファイルに作成した:
clear; function [val lev] = savecontour(filename,cnt) fid=fopen(filename,'w'); npar = 0; ptr = 0; N = size(cnt,2); while(ptr < N) npar++; ptr++; lev{npar} = cnt(1,ptr); n = cnt(2,ptr); val{npar} = cnt(:,((ptr+1):(ptr+n)))'; for i=1:n fprintf(fid,"%f %f\n",cnt(1,ptr+i),cnt(2,ptr+i)); end; fprintf(fid,"\n"); ptr+=n; endwhile; fclose(fid); endfunction; #------------- A=load("../jacobi_100.dat");//あんたがC++で作成したデータファイルを指定 nx=101; ny=101; x=A(:,1); y=A(:,2); u=A(:,3); x=reshape(x,ny,nx); y=reshape(y,ny,nx); u=reshape(u,ny,nx); [ cnt lev ]=contourc(x,y,u,7); [ val lev ]=savecontour("contour_100.dat",cnt);//等高線をcontour_100.datに保存
作成したコードは, Octave画面で次のように実行できる:
>> run jacobi.m
これで contour_100.dat ができたので, 作図した:
graph -Tps -C -x 0 1 -y 0 1 -nx -ny -X '\fIx' -Y '\fIy'\ --label-font-name Times-Italic --font-size 0.055\ --frame-line-width 0.0025 --line-width 0.0025\ contour_100.dat |ps2pdf -dEPSCrop - contour_100.pdf
これで contour_100.pdf が完成する.ここまで graph コマンドが長いと, draw.sh の利用を検討した方がいいかもしれない.