B. Mathematica に持っていって図を描く

MATLABで計算したデータを Mathematica に持っていって、 そちらで図を描いてみる。

• MATLAB で save('ファイル名', '変数名', '-ascii') として保存し た2次元配列は、 (1行の文字数が多くなり、 とても普通のテキスト・ファイルとは言えない代物だけれど) Mathematica では a=Import["ファイル名"]; で読める。
• MATLABでは長方形の各辺を $1280$ 等分して作った格子点上の 数値データが普通に扱えるが、 Mathematica には「想定外」の重さとなってしまうようだ (それが原因か定かでないが、Mathematica が落ちたことがある)。 計算その物は細かい格子を使うにしても、 描画の際にはデータを間引くのが良いかもしれない。
  $$N=10,20,40,80,160,320,640,1280,\dots$$
  なので、当面 $160$ 等分にしてみる。

  N=320
  [v320,d320]=eigs(plate_f1(N,0.3),200,0);
  

  あるいは

  load('free_0.3_320.mat')
  

  第4固有関数を抜き出して変換する。

  u=nv_to_dim2(v320(:,201-4),N,N);
  

  この u を間引いて $161\times161$ のサイズにする。

  i=1:N/160:N+1;
  j=1:N/160:N+1;
  uu=u(i,j);
  

  x=0:1/160:1;
  y=0:1/160:1;
  clf
  contour(x,y,uu)
  pbaspect([1 1 1])
  

間引き方の練習。

% drawcont.m
%  clear
%  load('free_0.3_1280.mat')
%  drawcont(v1280(:,201-n));figure(1)
function drawcont(v)
  [N,dummy]=size(v);
  N=sqrt(N)-1;
  if mod(N,160)==0
      u=nv_to_dim2(v,N,N);
      i=1:N/160:N+1;
      j=1:N/160:N+1;
      uu=u(i,j);
      x=0:1/160:1;
      y=0:1/160:1;
      clf
      contour(x,y,uu)
      pbaspect([1 1 1])
  else
     'N is not a multuple of 160'
  end

Mathematica で読むことを前提に一つ一つばらす。

clear
load('free_0.3_1280.mat')
i=1:N/160:N+1;
j=1:N/160:N+1;
for n=1:200
  u=nv_to_dim2(v1280(:,201-n),1280,1280);
  u=u(i,j);
  save(['poisson0.3/u' int2str(n) '.dat'], 'u', '-ascii')
end
e=diag(d1280);
e=e(200:-1:1);
save('poisson0.3/eigen.dat','e','-ascii')

これで poisson0.3 に 'u*.dat' と eigen.dat が 得られる。

u7=Import["Documents/MATLAB/poisson0.3/u7.dat"];
u8=Import["Documents/MATLAB/poisson0.3/u8.dat"];
ListPlot3D[u7,DataRange->{{0,1},{0,1}}]
ListPlot3D[u8,DataRange->{{0,1},{0,1}}]
ListContourPlot[u7]
ListContourPlot[u8]

Manipulate[ListContourPlot[a*u7+b*u8,DataRange->{{0,1},{0,1}},
                                   Contours->{0},ContourShading->None],
   {{a,1,"a"},-1,1,0.001},{{b,0.658,"b"},-1,1,0.001}]

For[i = 1, i <= 200, i++, 
 u[i] = Import[
   "Documents/MATLAB/poisson0.3/u" <> IntegerString[i] <> ".dat"]]

blend[m_, n_] := 
 Manipulate[
  ListContourPlot[a*u[m] + b*u[n], DataRange -> {{0, 1}, {0, 1}}, 
   Contours -> {0}, ContourShading -> None], {{a, 1, "a"}, -1, 1, 
   0.001}, {{b, 0.658, "b"}, -1, 1, 0.001}]

blend2[m_, n_] := 
 Manipulate[
  ListContourPlot[a*u[m] + b*u[n], DataRange -> {{0, 1}, {0, 1}}],
   {{a, 1, "a"}, -1, 1, 0.001}, {{b, 0.658, "b"}, -1, 1, 0.001}]

blend3[m_, n_] := 
 Manipulate[{ListPlot3D[a*u[m] + b*u[n], 
    DataRange -> {{0, 1}, {0, 1}}], 
   ListContourPlot[a*u[m] + b*u[n]]}, {{a, 1, "a"}, -1, 1, 
   0.001}, {{b, 0.658, "b"}, -1, 1, 0.001}]