XLPack for Python および XLPack for Matplotlib の使用例 (1) 3D 表示 (ラプラス方程式)

\(x = 0 \sim 1, y = 0 \sim 1\) の正方形領域において, 次の偏微分方程式を解くものとします.
\[
-d^2u/dx^2 – d^2u/dy^2 = 0
\] ただし, 境界条件は
\[
u(x, 0) = 0, u(0, y) = 0, u(x, 1) = x, u(1, y) = y \space (ディリクレ条件)
\] とします.

このような偏微分方程式は, 例えば熱平衡状態にある板の温度分布の記述などに適用されます.

この方程式の解法としては差分法がよく使用されます. 5 点差分近似により数値解を求める例は以下を参照してください.

16. 疎行列の線形計算

ここでは有限要素法を用いて解を求め 3D 表示する例を説明します.

\(x, y\) 方向に正方形領域をそれぞれ \(Nx, Ny\) 等分して長方形格子に分割します. さらにそれぞれを２分割して三角要素とした単純メッシュを生成することにします.

これを用いて有限要素法により解を求めます. なお, 上図は 7 等分 (格子点数 8 x 8)の例を示しますが, 実際の計算では 19 等分 (格子点数 20 x 20) (Nx = 19, Ny = 19) としました.

1. 全体の構成

メインプログラムは以下のとおりです.

Sub Main()
    Const Nx = 19, Ny = 19
    Const Nb2 = 0, LdKnc = 4, LdKs = 3
    Const N = (Nx + 1) * (Ny + 1), Ne = 2 * Nx * Ny, Nb = 2 * (Nx + Ny)
    Dim X(N - 1) As Double, Y(N - 1) As Double
    Dim P(N - 1) As Double, Q(N - 1) As Double, F(N - 1) As Double
    Dim Knc(LdKnc - 1, Ne - 1) As Long, Ks(LdKs - 1, Nb - 1) As Long
    Dim Lb(Nb - 1) As Long, Ib(Nb - 1) As Long, Bv(Nb - 1) As Double
    Dim Ks2() As Long, Alpha() As Double, Beta() As Double
    Dim Val(20 * N) As Double, Ptr(N) As Long, Ind(20 * N) As Long
    Dim B(N - 1) As Double, U(N - 1) As Double
    Dim Ux As Double, Err As Double
    Dim Info As Long
    Dim I As Long, J As Long
    '-- Set mesh data
    Call SetData(Nx, Ny, N, Ne, Nb, X(), Y(), Knc(), Ks(), Lb(), P(), Q(), F(), Ib(), Bv())
    '-- Assemble FEM matrix
    Call Fem2p(N, Ne, X(), Y(), Knc(), P(), Q(), F(), Nb, Ib(), Bv(), Nb2, Ks2(), Alpha(), Beta(), Val(), Ptr(), Ind(), B(), Info)
    '-- Solve equation by CG method
    Dim Iter As Long, Res As Double
    Call Cg1(N, Val(), Ptr(), Ind(), B(), U(), Info, Iter, Res)
    Debug.Print "Info =" + Str(Info) + ", Iter =" + Str(Iter) + ", Res =" + Str(Res)
    If Info <> 0 Then Exit Sub
    '-- Plot solution
    Call Plot(Nx, Ny, U())
End Sub

XLPack ライブラリのサブルーチン Fem2p() は有限要素行列の組み立てを行います. 疎行列 (CSR) 形式の係数行列と右辺ベクトルからなる方程式が得られます. 係数行列は対称行列なので, 方程式は XLPack ライブラリのサブルーチン Cg1() を使って共役勾配法により解くことができます.

サブルーチン SetData() は Fem2p() の入力となる格子点, 有限要素, および, 境界条件のデータを設定します.

得られた解はサブルーチン Plot() により表示されます. Plot() では XLPack for Python または XLPack for Matplotlib を使用して 3D 表示を行います.

実行結果は次のようになります.

2. 入力データの設定

三角要素による単純メッシュは XLPack ライブラリのサブルーチン Mesh23() を使って生成することができます.

次に, ラプラス方程式の係数 p(x, y) = 1, q(x, y) = 0, f(x, y) = 0 を設定します. また, ディリクレ境界条件を設定します.

これらを行うのが次のサブルーチンです.

Sub SetData(Nx As Long, Ny As Long, N As Long, Ne As Long, Nb As Long, X() As Double, Y() As Double, Knc() As Long, Ks() As Long, Lb() As Long, P() As Double, Q() As Double, F() As Double, Ib() As Long, Bv() As Double)
    Dim Id() As Long
    Dim I As Long, J As Long, K As Long
    '-- Generates simple rectangular mesh for 2D triangular element
    Call Mesh23(Nx, Ny, X(), Y(), Knc(), Ks(), Lb())
    '-- Set p(x, y) = 1, q(x, y) = 0 and f(x, y) = 0
    For I = 0 To N - 1
        P(I) = 1
        Q(I) = 0
        F(I) = 0
    Next
    '-- Set boundary condition 1
    ReDim Id(N - 1)
    For I = 0 To Nb - 1
        For J = 1 To 2
            Id(Ks(J, I) - 1) = Lb(I)
        Next
    Next
    K = 0
    For I = 0 To N - 1
        If Id(I) = 1 Or Id(I) = 4 Then
            Bv(K) = 0
        ElseIf Id(I) = 2 Then
            Bv(K) = Y(I)
        ElseIf Id(I) = 3 Then
            Bv(K) = X(I)
        Else
            GoTo Continue
        End If
        Ib(K) = I + 1
        K = K + 1
Continue:
    Next
End Sub

3. 解の 3D 表示

得られた解を 表示する Python のコードを XLPack for Python を使って直接記述することができます.

Sub Plot(Nx As Long, Ny As Long, U() As Double)
    Dim I As Long
    For I = 0 To (Nx + 1) * (Ny + 1) - 1
        Cells(6 + I, 2) = U(I)
    Next
    '-- Python codes
    PY "import matplotlib.pyplot as plt"
    PY "from matplotlib import cm"
    PY "import numpy as np"
    PY "from XLPackPy import *"
    PY "nx =" + Str(Nx) + "+ 1"
    PY "ny =" + Str(Ny) + "+ 1"
    PY "U = xlc_get(5, 1, 4 + nx*ny)"
    PY "X = np.empty((ny, nx))"
    PY "Y = np.empty((ny, nx))"
    PY "Z = np.empty((ny, nx))"
    PY _
      "for i in range(ny):" + vbNewLine + _
      "  for j in range(nx):" + vbNewLine + _
      "    X[i][j] = (1.0/(nx - 1))*j" + vbNewLine + _
      "    Y[i][j] = (1.0/(ny - 1))*i" + vbNewLine + _
      "    Z[i][j] = U[nx*i + j]"
    PY "fig = plt.figure()"
    PY "ax3 = fig.add_subplot(projection='3d')"
    PY "surf = ax3.plot_surface(X, Y, Z, cmap=cm.coolwarm)"
    PY "fig.colorbar(surf, shrink=0.5, aspect=5)"
    PY "plt.show()"
End Sub

XLPack for Matplotlib を使うと全て VBA で記述することができます. Fig.Add_subplot_3d() は Python の fig.add_subplot(projection=’3d’) に相当します. その他は Python のコードとほぼ同じになります.

Dim Plt As New Pyplot

Sub Plot(Nx As Long, Ny As Long, U() As Double)
    Dim X() As Double, Y() As Double, Z() As Double
    ReDim X(Ny, Nx), Y(Ny, Nx), Z(Ny, Nx)
    Dim I As Long, J As Long
    Dim Fig As Figure, Ax3 As Axs3d, Surf As PyObject
    '-- Prepare data
    For I = 0 To Ny
        For J = 0 To Nx
            X(I, J) = (1# / Nx) * I
            Y(I, J) = (1# / Ny) * J
            Z(I, J) = U((Nx + 1) * I + J)
        Next
    Next
    '-- Plot by Matplotlib
    Set Fig = Plt.Figure()
    Set Ax3 = Fig.Add_subplot_3d()
    Set Surf = Ax3.Plot_surface(Ny + 1, Nx + 1, X(), Y(), Z(), "coolwarm")
    Call Fig.Colorbar(Surf, KwArgs:="shrink=0.5, aspect=5")
    Call Plt.Show
End Sub
公開日: 2026/04/28 (更新日: 2026/04/28)