howto/maxima †

Maximaをベースにした種類 †

wxmaxima, 12.04.0(win7)
├ maxima, 5.27.0(win7)
└ xmaxima, テキストベースらしい

imaxima, wxmaximで十分のようだ

texmacs, 1.99.1.8496(win7), maximaとの連携がうまくいかない。それ以外は単なるエディタ。

• mac版のwxmaximaはファイル保存できない。.wxmxと.macには入力テキストいちおう保存される(2013年)

一行メモ †

• オイラーの公式を展開する例（これしかない？

  trigsimp(
    subst( cos(phi)+%i*sin(phi), exp(%i*phi), 
     1/sqrt(2)/%i*( spherical_harmonic(1,1,theta,phi) +  spherical_harmonic(1,-1,theta,phi) )
    )
  );

• 水素原子の波動関数

  LaguerreBiPoly(k,s,x):= diff( exp(x)*diff(x^k*exp(-x),x,k), x, s );
  LaguerreBiFunc(n,l,x):= expand( exp(-x/2)*x^l * LaguerreBiPoly(n+l,2*l+1,x) );
  SpheHarm(l,m,theta,phi):=spherical_harmonic(l,m,theta,phi);   /* 2014/11/20訂正,2016/07/13訂正 */
                             /* このsperical_harmonic には、通常のものに(-1)^(m)が掛けられている */
  SpheHarm(l,m,theta,phi):=(-1)^(m)*spherical_harmonic(l,m,theta,phi); /* 通常のYlm */
  Psi(n,l,m,r,theta,phi):= -sqrt( (n-l-1)!/(2*n*((n+l)!)^3)*(2/n)^3 ) *
    factor(( subst( 2*r/n, x, LaguerreBiFunc(n,l,x) ) * SpheHarm(l,m,theta,phi) ));

• auからfsecへの変換

  fs2au(t):=t/2.418e-2;

• Maxima 5.33.0 マニュアル（日本語） http://maxima.sourceforge.jp/maxima.html -- 2014-06-05 (木) 10:15:22
• trgsimp以外に、radcanという指数・対数の整理があり。 -- 2014-06-05 (木) 10:49:54
• declare(k,integer) k を整数として宣言 -- 2014-06-05 (木) 10:54:23
• 前置演算子 -- 2014-06-05 (木) 10:56:00

  prefix( "A" );
  A f := diff( r^2 * diff( f, r ), r )/(2*r^2)
       + diff( sin(theta) * diff( f, theta ), theta )/(2*r^2*sin(theta))
       - 1/r*f - E*r*sin(omega*t)*cos(theta)*f;
  Psi[0](r):=exp(-r)/sqrt(%pi);
  A Psi[0](r);

• 部分分数分解 partfrac( (2*x^2+6*x+1)/(x*(x+1)^2), x ); -- 2014-06-19 (木) 12:14:31
• 式をまとめる？： factorout, rat, ordergreat -- 2014-09-10 (水) 10:57:54
• http://www.kabipan.com/computer/maxima/quick.html より -- 2014-09-14 (日) 20:21:46

  1/3, numer; と入力すると .3333333333333333 が得られる
  
  値の割り当て（assingment）
  シンボルはその場で評価され、その値に置き換わる
  remvalue (n); で、n に割り当てられた値を削除する。
  remvalue (all); で、ユーザがシンボルに割り当てた値をすべて削除する。
  kill(n); でシンボル n を削除。kill (all) で、ユーザが作成したすべてのシンボルを削除。
  シンボルの値として割り当てることができるものは、数に限らない。たとえば式を値として
  割り当てることができる：シンボルに割り当てられた値を評価したい場合は、ev () を用いる
  
  declare (n, integer); で、n を整数であると宣言する。
  facts (); で n 以外の変数に対する制約も表示される。
  フラグは features; で知ることができる
  remove (all,integer); で、すべてのシンボルから整数フラグを取り除くことができる
  
  assume (n > 0); で、n が 0 より大きいことを仮定することができる
  forget (n > 0) で、仮定 n > 0 を取り除くことができる。forget (all) で、すべての assume されたものを削除
  
  fundef(f); で関数 f の定義を見ることができる。dispfun(all) ですべてのユーザ関数定義を見ることができる。
  remove (f, function); で、f から関数定義を取り除くことができる
  
  assume (n = m); の形が許されないのも、これに似た事情である。n の字面はどうやっても m の
  字面とは同じではない。assume (equal(n, m)); なら可能である。

• 行列をくっつけるには、行の場合はaddcol(A, B)、列の場合は addrow(A, B) -- 2014-12-10 (水) 14:57:34
• 行列の要素を取り出すには、Eq[i,j]である。一列や一行であってもi,jを指定する。.jを略すとi行目の一行の行列が得られるようだ -- 2014-12-10 (水) 15:03:27
• x[k]の扱い方。（配列、配列関数とかいうものらしい）

  x[2*k] は 下付き添字2kが付いた x で表示されるが、配列のようだ。
  arrayinfo( x ) で、[hashed, 次元, 登録されたハッシュキー, …]
  listarray( x) で、[ハッシュ値, …]
  を知ることができる。kill( x ) も可。
  
  下付き文字 2k の付いた x の文字として使用するには、' を付ける。
  例） integrate( 'x[2*k+1], x, 'x[2*k+1]-h, 'x[2*k+1]+h );

• 配列関数 N[9](x,y):=1/1*(1-x)*(1+x)*(1-y)*(1+y); -- 2014-12-20 (土) 09:53:46
• h[i,j]:=(i-1)*3+j; genmatrix(h,3,3); で行列を作成できる。 -- 2014-12-20 (土) 16:43:45
• row(A,2); col(A,2); 行列の行や列を取り出す。 -- 2014-12-20 (土) 16:44:42
• setelmx(p, 2, 3, A); 行列Aの2行3列をpに置き換える。 -- 2014-12-20 (土) 16:45:40
• a: lambda([i, j], x[j]ˆi); genmatrix(a, 3, 4, 0, 1); ラムダを使った行列の作成 -- 2014-12-20 (土) 16:49:46
• 連立方程式を解く -- 2014-12-25 (木) 11:00:27

  B0(x):=a*x^5+b*x^4+c*x^3+d*x^2+e*x+f;
  B1(x):=a*5*x^4+b*4*x^3+c*3*x^2+d*2*x+e;
  B2(x):=a*5*4*x^3+b*4*3*x^2+c*3*2*x+d*2;
  /* [-1,0] j=0 */
  solve( [B0(0)=1,B1(0)=0,B2(0)=0,B0(-1)=0,B1(-1)=0,B2(-1)=0],[a,b,c,d,e,f] );

• 極限値 limit(sin(x)/x, x, inf); limit(1/x, x, 0, minus); 4番目の引数は方向を示す -- 2015-02-10 (火) 17:47:03
• 等高線を描く wxcontour_plot ( uo(x)*uo(y), [x,-1,1], [y,-1,1] ); -- 2015-02-28 (土) 10:42:25
• 数値積分 romberg( exp(-x^2), x, 0, 1 ); rombergtol;(なおこれはrombergを使用したら値が初期値がセットされる); 被積分関数に条件分岐(if)があればやはり積分してくれない。 -- 2015-03-03 (火) 18:04:22
• 積分の例： integrate( r*r*exp(-2*r) , r, 0, inf ); -- 2016-07-05 (火) 10:47:54

ベクトルの外積 †

関数を定義して使うしかないらしい
cross( u, v ):=  [ u[2]*v[3] - u[3]*v[2], u[3]*v[1] - u[1]*v[3], u[1]*v[2] - u[2]*v[1] ];
A1:[a,0,0] ;
A2:[0,b,0];
A3:[0,0,c];
A1 . cross( A2, A3 ); /* スカラー三重積 */

グラフ描画 †

• plot2d

  f(x, A, a):=x*(1+A*exp(-a*x));
  plot2d( [f(x,4,2), f(x,8,2),f(x,4,4),f(x,8,4), x], [x, 0,5],
    [legend, "f(x,4,2)", "f(x,8,2)","f(x,4,4)", "f(x,8,4)", "x"],
    [gnuplot_preamble, "set grid; set key bottom right;
                        set xrange[0:5]; set yrange[0:5];"] );
  またはwxplot2d
  
  軸ラベルは [xlabel,"Intensity (W/m^2)"], [ylabel,"E (V/Angstrom)"]
  軸を対数にするには、[gnuplot_preamble, " set logscale x; set logscale y; "]
  軸の目盛の特殊なフォーマットは、[gnuplot_preamble, "set format x \"%1.0le%2L\"; "] Lは10のべき乗部分、lは仮数部
  
  別の例)
  wxplot2d( [D(t/au2fs),C1(t/au2fs),C2(t/au2fs)], [t,0,1086*au2fs],
      [legend, "Dipole", "C1", "C2"],
      [xlabel,"time (fs)"], [ylabel,"(au)"],
      [gnuplot_preamble, "set grid; set key top center;" ]
  );

• plot3d

  f(x,t):=-400/16*sin(2*%pi/0.4*t)*(x-8);
  wxplot3d( f(x,t), [x,0,16], [t,0,0.4],
      [legend, "f(x,t)"] );
  またはwxplot3d

• plot2d（媒介変数）

  plot2d( [parametric, 1/2*t*t+t, sin(2*%pi*t)], [t,0,50], [nticks, 10000] );
  変数はtでないとダメなようだ。

• その他

  [legend,false]  凡例を表示させない

雑多なこと †

assume( W > 0 )  /* 前提条件。積分とかで場合分けの質問が煩わしいときなどに使用。andを使用可。 */
forget( W )  /* 変数Wの値を未設定にする？ */
forget( all )  /* 全ての変数の値を未定義にする？ */
kill( all )  /* 全ての変数の値を未定義にする？ */

W:0.25         /* 変数Wの値を0.25に設定する */

inf  /* ∞ */
%i   /* 虚数単位 */
%pi  /* π */

conjugate( hoge )  /* 複素共役をとる */
realpart( hoge )   /* 実部を取り出す */
imagpart( hoge )   /* 虚部を取り出す*/
abs( hoge )        /* 絶対値 */

diff( hoge(x), x [,n] )  /* hoge(x)のn回微分。nを省略したときは1 */

taylor( exp(%i*C*x), x, 0,5)  /* テイラー展開をxについて、０次から５時までする */

subst( a, x, hoge(x) )    /* 関数hogeのxをaに置き換える（x,aは数式でも可） */

nusum( A^^k, k, 0, inf );
num( BBconjinv[1,1] );

末尾$    /* 出力を省略 */

hoge(x), x:0.25;  /* env( x:0.25, hoge(x) ) と同じく、一時的な代入 */

romberg  数値積分

関数の定義 †

hoge( x, y, z ):= x,y,zを使った関数  /* 関数の定義 */

f(x):=0.3*x*x-1.2*x+5;
define( nab_f(x), diff( f(x), x ) );  /* 微分として定義する場合 */

nab_f(x):=diff( f(x), x );    /* 期待通りに働かない */
nab_f(x):='diff( f(x), x );   /* 期待通りに働かない */

f(1.2);
nab_f(1.2);

場合分けのある関数の定義
ua(x,c):=if x>=c-1 and x<=c+1 then 1/2*(x-c)*(x-c-1) else 0;

行列 †

A:matrix( [ 1,2],[3,4] );
A:matrix( [-1, exp(-%i*phi)], [exp(%i*phi), 0] );
A.A; A^^3; A^^4; A^^5; A^^-1; A^^0;  /*  行列の積は . べき乗は ^^ 逆行列は^^-1 */

exp( A ); /* これは各要素のexpが計算されるので注意 */

determinant(A);  /* 行列式の値。なお、0のときA^^-1はエラーになる */

transpose( A )  /* 行列Aを転置する */
conjugate(transpose( A ))  /* Aの転置複素共役 */

浮動点少数として計算 †

float( 1/3 );  /* 浮動点少数として結果を出力 */
bfloat(1/3);   /* 〃 */
bfloat(1/3), fpprec:8;  /* 8桁の浮動点少数として出力 */

方程式を解く †

TOUfunc(V0,E,a):=( (1+V0^2*(sinh(sqrt(2*(V0-E))*a))^2/(4*E*(V0-E)) )^(-1) );
ans:solve( TOUfunc( 80, 74, a )=0.5, a );
bfloat( [rhs(ans[1]), rhs(ans[2])] ), fpprec:4;

式展開やまとめる方法 †

expand   展開 <--> factor 因数分解
partfrac      <--> ratsimp
trigreduce 三角関数の積を減少   <--> trigsimp
trigexpand 三角関数の引数を展開 <--> trigsimp