数式処理

Matlabのsymbolic math toolboxの備忘録です。

表記法
sym(4/3,'r');%有理数
4/3
sym(pi,'r') %有理数
pi
sym(pi,'f'); %浮動小数点 1.F*2^E; F 16進数 E 指数
1.921fb54442d18'*2^(1)
digits(10);%少数第10桁まで
sym(pi,'d');% 少数表示
3.141592654
vpa(pi,20); % 20桁表示
3.1415926535897932385
sym(1e-10,'d');
.1000000000e-9 % 1.0000000000e-10にして欲しいなあ
sym('1+sqrt(3)'); % √表示
1+sqrt(3) % 文字列のままです
simplify(sym('1+sqrt(3)'))
1+3^(1/2) % √表示
x=sym('(1+2*sqrt(3))*(2-3*sqrt(3))');
simplify(x) % √表示
-(1+2*3^(1/2))*(-2+3*3^(1/2))
expand(x) % 展開
x=
-16+3^(1/2)
simplify(sym('1+exp(2)+sin(3)+cos(pi)'))
exp(2)+sin(3) % 略せない関数はそのまま残されます

simplify(sym(文字列))として書くようにしたほうがよさそうです。

LaTeX表示
syms a b c x;
s=solve('a*x^2+b*x+c=0',x)%ax^2+bx+c=0を解く
s=
-1/2*(b-(b^2-4*a*c)^(1/2))/a
-1/2*(b+(b^2-4*a*c)^(1/2))/a
latex(s) %LaTeXに変換　1行で表示されるのですが、長いので適当に改行しています。
\left[ \begin {array}{c}
-1/2\,{\frac {b-\sqrt {{b}^{2}-4\,ac}}{a}}\\
\noalign{\medskip}
-1/2\,{\frac {b+\sqrt {{b}^{2}-4\,ac}}{a}}
\end {array} \right]
Matlabのコマンドウインドウ上にはLaTeXは反映されません。figure上で反映させます。
LaTeXの数式を表す$・・・$の形式にするため strcat('$',latex(s),'$')という構文にします。 \[や\begin{equation}などは使えません。 figure;text(0,0.5,strcat('$',latex(s),'$'),'interpreter','latex','fontsize',18);

プログラムのしやすさを考慮し、見た目の分かりやすさは考慮していないLaTeX文になっています。
LaTeXの文法に基づいて以下のように手動で書き換えることもできます。
ss='\left[ \begin{array}{c} -\frac{b-\sqrt{b^2-4ac}}{2a} \\ -\frac{b+\sqrt{b^2+4ac}}{2a} \end{array} \right]'
figure;text(0,0.5,strcat('$',ss,'$'),'interpreter','latex','fontsize',18);


代入
syms x y;f=x^2+y^2;
subs(f,[x,y],[1,3]) % x=1 y=3を代入
10
syms x y t real; % 実数のみとする
f=cos(x)^2+sin(y)^2;
k=subs(f,[x,y],[t,t]) % x=t y=tを代入
k=
cos(t)^2+sin(t)^2
simplify(k) %簡略化
1
expand(k) %展開
cos(t)^2+sin(t)^2 % 簡略化はしてくれません


連立方程式

syms X Y Z DetA A11 A12 A13 A21 A22 A23 A31 A32 A33 x y z;
Eq1='X=A11*x+A12*y+A13*z'; %文字列表記の方が後々楽です
Eq2='Y=A21*x+A22*y+A23*z';
Eq3='Z=A31*x+A32*y+A33*z';
[x,y,z]=solve(Eq1,Eq2,Eq3)
x=
(-X*A32*A23+X*A22*A33-A12*Y*A33+A12*A23*Z-A13*A22*Z+A13*A32*Y)/(-A11*A32*A23+A11*A22*A33-A22*A31*A13-A21*A12*A33+A32*A21*A13+A31*A12*A23)
y=
-(-A11*Y*A33+A11*A23*Z-A21*A13*Z-A23*A31*X+Y*A31*A13+A21*X*A33)/(-A11*A32*A23+A11*A22*A33-A22*A31*A13-A21*A12*A33+A32*A21*A13+A31*A12*A23)
z=
(A11*A22*Z+A31*A12*Y-A22*A31*X-A11*A32*Y-A21*A12*Z+A32*A21*X)/(-A11*A32*A23+A11*A22*A33-A22*A31*A13-A21*A12*A33+A32*A21*A13+A31*A12*A23)
[r,s]=subexpr([x;y;z],'s') % 共通項sを見つけ、sでまとめる
r =
conj((A12*A23*Z-A12*A33*Y-A13*A22*Z-A23*A32*X+A33*A22*X+A13*A32*Y)/s)
-conj((X*A21*A33-X*A31*A23+A11*A23*Z-A11*A33*Y-A13*A21*Z+A13*Y*A31)/s)
conj((A22*A11*Z-A22*X*A31+X*A21*A32-A11*A32*Y-A12*A21*Z+A12*Y*A31)/s)
s=
A11*A22*A33-A11*A23*A32-A21*A12*A33+A21*A13*A32+A31*A12*A23-A31*A13*A22


連立方程式を配列表示（非推奨）

syms Q A P;
Eq1='Q(1)=A(1,1)*P(1)+A(1,2)*P(2)+A(1,3)*P(3)';
Eq2='Q(2)=A(2,1)*P(1)+A(2,2)*P(2)+A(2,3)*P(3)';
Eq3='Q(3)=A(3,1)*P(1)+A(3,2)*P(2)+A(3,3)*P(3)';
[P(1),P(2),P(3)]=solve(Eq1,Eq2,Eq3);
警告: 3 equations in 1 variables.
> In solve at 113
警告: Explicit solution could not be found.
> In solve at 140
??? 代入文 A(:) = B では、A と B の要素数が同じである必要があります
求める変数を配列にすると解けないみたい。

syms Q A x y z;
Eq1='Q(1)=A(1,1)*x+A(1,2)*y+A(1,3)*z';
Eq2='Q(2)=A(2,1)*x+A(2,2)*y+A(2,3)*z';
Eq3='Q(3)=A(3,1)*x+A(3,2)*y+A(3,3)*z';
[x,y,z]=solve(Eq1,Eq2,Eq3)
x=
(A(1,2)*Q(3)*A(2,3)-A(1,2)*A(3,3)*Q(2)-A(1,3)*A(2,2)*Q(3)+A(1,3)*Q(2)*A(3,2)+Q(1)*A(2,2)*A(3,3)-Q(1)*A(2,3)*A(3,2))/(-A(2,1)*A(1,2)*A(3,3)-A(2,2)*A(3,1)*A(1,3)+A(2,2)*A(1,1)*A(3,3)+A(2,1)*A(1,3)*A(3,2)+A(2,3)*A(3,1)*A(1,2)-A(2,3)*A(1,1)*A(3,2))
y=
-(A(1,1)*Q(3)*A(2,3)-A(1,1)*A(3,3)*Q(2)+A(3,1)*A(1,3)*Q(2)-A(3,1)*Q(1)*A(2,3)+A(3,3)*A(2,1)*Q(1)-Q(3)*A(2,1)*A(1,3))/(-A(2,1)*A(1,2)*A(3,3)-A(2,2)*A(3,1)*A(1,3)+A(2,2)*A(1,1)*A(3,3)+A(2,1)*A(1,3)*A(3,2)+A(2,3)*A(3,1)*A(1,2)-A(2,3)*A(1,1)*A(3,2))
z=
1/(-A(2,1)*A(1,2)*A(3,3)-A(2,2)*A(3,1)*A(1,3)+A(2,2)*A(1,1)*A(3,3)+A(2,1)*A(1,3)*A(3,2)+A(2,3)*A(3,1)*A(1,2)-A(2,3)*A(1,1)*A(3,2))*(-A(2,1)*A(1,2)*Q(3)+A(2,2)*A(1,1)*Q(3)-A(2,2)*A(3,1)*Q(1)+Q(2)*A(3,1)*A(1,2)-Q(2)*A(1,1)*A(3,2)+A(2,1)*Q(1)*A(3,2))

Mathematicaと異なり、基本的に記号symbolicの配列定義はちょっと面倒です。
A=repmat(sym(0),3,3)
と予め定義しておくほうがよさそうです。

Biot-Savar
syms Qx Qy Qz ax ay az t;
B=(4*pi*1e-7)/(4*pi)*cross([Qx;Qy;Qz],[ax;ay;az])/(ax^2+ay^2+az^2)^(3/2)*1e+15;
B =simplify(B)
100000000*(Qy*az-Qz*ay)/(ax^2+ay^2+az^2)^(3/2)
100000000*(Qz*ax-Qx*az)/(ax^2+ay^2+az^2)^(3/2)
-100000000*(-Qx*ay+Qy*ax)/(ax^2+ay^2+az^2)^(3/2)
Bz=subs(B(1),[Qx,Qy,Qz,ax,ay,az],[0,0,1e-9,0.05,0.01*t,0]);% [0,0,1]nA, [5,x,0]cm上
figure('color',[1,1,1]);ezplot(Bz,[-20,20]);grid on;

Bz=subs(B(1),[Qx,Qy,Qz,ax],[0,0,1e-9,0.05]);
ezsurf(Bz,[-1,1]*0.5);