MATLAB,Python,Scilab,Julia比較 第2章 その36【対称行列と二次形式⑧】

バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia2-backnumber/

はじめに
登場人物
ロードマップおよび数式【再掲】
Juliaコード
実行結果
考察
まとめ

はじめに

正規方程式を導出するまでの説明。
今回は二次形式の多項式表現と行列表現の計算をJuliaで実施する。

登場人物

博識フクロウのフクさん

イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1

エンジニア歴8年の太郎くん

イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1

ロードマップおよび数式【再掲】

太郎くん

ロードマップと今回使用する数式を再掲。

正規方程式に至る道、二次形式、対称行列、二次形式の微分、グラム行列、二乗和誤差、正規方程式

\(
x^TA\circ x^T
\begin{bmatrix}
1\\
\vdots\\
1
\end{bmatrix}
\)

計算する二次形式の数式

\(
3x^2+2y^2+5xy
\)

フクさん

今回は、Juliaで実施する。

Juliaコード

フクさん

Juliaコードは以下になる。

using PyPlot

function meshgrid(xin,yin)
	nx=length(xin)
	ny=length(yin)
	xout=zeros(ny,nx)
	yout=zeros(ny,nx)
	for jx=1:nx
	    for ix=1:ny
	        xout[ix,jx]=xin[jx]
	        yout[ix,jx]=yin[ix]
	    end
	end
	return (x=xout, y=yout)
end

a=3;
b=2;
c=5;
A=[a c/2; c/2 b];

N=6;
ax=range(0,10,length=N);
ay=range(0,10,length=N);
x,y=meshgrid(ax,ay);

fig, (ax1, ax2) = plt.subplots(2, 1, 
	figsize=(8, 8), 
	subplot_kw=Dict("projection" => "3d"))

polyY=a.*x.^2+b.*y.^2+c.*x.*y;
ax1.plot_wireframe(x, y, polyY, rstride=1, cstride=1);

X=[x[:]';y[:]'];
Y=X'*A.*X'*ones(2,1);
matY=reshape(Y,N,N);

ax2.plot_wireframe(x, y, matY);

plt.tight_layout()
plt.show()

polyY
matY

実行結果

フクさん

実行結果は以下になる。

polyY
6×6 Matrix{Float64}:
   0.0   12.0   48.0  108.0  192.0   300.0
   8.0   40.0   96.0  176.0  280.0   408.0
  32.0   84.0  160.0  260.0  384.0   532.0
  72.0  144.0  240.0  360.0  504.0   672.0
 128.0  220.0  336.0  476.0  640.0   828.0
 200.0  312.0  448.0  608.0  792.0  1000.0

matY
6×6 Matrix{Float64}:
   0.0   12.0   48.0  108.0  192.0   300.0
   8.0   40.0   96.0  176.0  280.0   408.0
  32.0   84.0  160.0  260.0  384.0   532.0
  72.0  144.0  240.0  360.0  504.0   672.0
 128.0  220.0  336.0  476.0  640.0   828.0
 200.0  312.0  448.0  608.0  792.0  1000.0