面向對象：MATLAB的自定義類 [MATLAB]

　　這幾天剛剛開始學習MATLAB的面向對象編程。以前做的事情都是用MATLAB寫一些簡單的腳本或者函數，這方面MATLAB成熟的函數和直截了當的矩陣運算方法和語法都很容易上手，方便人專注於算法本身。前幾天寫代碼的時候想到，在實際用MATLAB進行數值計算時，將數據和函數用一些方法組織起來也會帶來很大的便利，否則零散的數據和函數總歸看着不舒服。比如，我恰好最近想在MATLAB里面寫一點代碼讓應力張量相關的計算變得簡潔一些。

一、應力張量示例的物理背景

　　在彈性力學里，應力張量是描述彈性體內部某一點的應力狀況的；它和過該點且法向量為 $\hat{n}$ 的面上的應力（矢量）關系是：$$\boldsymbol{t}=\boldsymbol{T}\cdot \hat{n}=\hat{n}\cdot \boldsymbol{T}$$　　這里粗寫的 $\boldsymbol{t}$ 是應力矢量；粗寫的 $\boldsymbol{T}$ 是應力張量，是一個二階張量，可以表示為一個三維矩陣。應力張量顯然不是一個可以直接觀測/測量的量，對於一般理解來說，應力矢量顯然是一個更為直觀、物理意義更加明確的物理量，因此已知應力張量求應力矢量應當是一個很基本的操作。

　　已知應力張量，求某一個面上的正應力和剪應力大小的方法是：$$\sigma_n=\hat{n}\cdot \boldsymbol{T}\cdot \hat{n},\quad \tau_n=\sqrt{|\boldsymbol{t}|^2-\sigma_n^2}$$　　正應力就是應力矢量和面的法向量再做一次內積，也可以看成是應力張量和法向量的二次內積。剪應力和正應力方向正交，其矢量和為應力矢量，因此用勾股定理就能得到其大小。在彈性材料中，剪應力往往是非常重要的一個物理量，因為許多材料（比如岩石等）都有抗壓不抗剪的特點，剪應力的大小將直接決定這些材料是否會發生破裂、如何（沿哪個方向）發生破裂。所以已知應力張量求正應力和剪應力也是很基本的一個操作。

　　此外，彈性體中的應力張量總是可以改寫成主軸坐標系的形式，其物理意義是總是存在三個正交的主應力方向，以這三個主應力方向為坐標軸表示，則應力張量剪切分量均為零；其數學意義是應力張量總是可以對角化。在主應力坐標軸下表示有特別的簡潔性，從而有一些好處，因此求出應力張量的主應力方向和主應力大小都具有重要的意義。

　　現在，我有了三個想要實現的操作。當然，因為MATLAB直觀人性化的矩陣運算語法，這三個操作在腳本里都很容易實現。比如，我現在有一個應力張量 T，三個操作分別可以用以下的代碼實現：

% 求應力矢量，要求為行向量
t = (T * n).';    % 若n為列向量
t = n * T;    % 若n為行向量

% 求正應力
sigma = n.' * T * n;    % 若n為列向量
sigma = n * T * n.';    % 若n為行向量
% 求剪應力
tau = sqrt(norm(t)^2 - sigma^2);

% 求主應力及其對應方向
[V, D] = eig(T);

　　看起來似乎也不復雜，但是這個存在四個不足之處：（一） 求應力矢量和求正應力時必須要判斷輸入的n是行向量還是列向量，如果每次寫一個判斷分支語句非常麻煩；（二）求剪應力必須以求正應力為基礎，這意味着每次求剪應力都需要手寫一遍求正應力的操作，這也包括判斷語句；（三）主應力的求解結果為兩個矩陣，如果想要直接得到“主應力-主應力方向”的結構怎么辦呢？（四）表達不直觀，可讀性比較差。

　　這些不足之處可以用自定義函數來解決，寫一個function把代碼放進去就可以了。但是這樣又有一個問題：這幾個操作還是沒有發生關聯！我們不知道這幾個操作都是對於同一個物理量——應力張量進行的操作。如果能像C/C++/Python那樣放個類來存放這些函數，然后再用T.X()這樣的方法調用豈不是美滋滋？這就是要學習——

 

 

二、定義應力張量的MATLAB類

　　現在，利用MATLAB中的“類”，將應力張量的幾個操作集合到一起。那么首先，需要定義一個MATLAB類：

% stressTensor.m
% 關鍵詞classdef，后跟一個自定義類名，存放屬性和方法，以end結束
classdef stressTensor
    %    應力張量(stressTensor)類

    % 關鍵詞properties：用來存放類的屬性（成員變量），以end結束
    properties
        tensorMat
    end

    % 關鍵詞methods：用來存放類的方法（成員函數），以end結束
    methods
        % 關鍵詞function：這個見過很多遍了，就是用來定義函數的，寫法和直接定義函數一樣，以end結束
        % 下為構造函數，名稱和類名稱相同，返回值必須是stressTensor類，因此在書寫時可以直接寫"item.tensorMat"，即訪問其作為stressTensor類具有的成員變量
        function item = stressTensor( A )
            if nargin == 0
                item.tensorMat = zeros(3, 3);
            else
                item.tensorMat = A;
            end
        end
    end
end

　　在定義一個類的時候，至少要定義一個構造函數（正如上文已經實現的），否則可能根本無法生成一個實例。這里有個小技巧，在函數中nargin是一個特殊的關鍵詞，它的涵義可以理解為：n-arguments-input，即傳入參數的個數。利用nargin，可以在函數結構內部判斷輸入了幾個參數，並且分別做出響應，這等於寫一個函數就達到了函數重載的作用。除此之外，MATLAB的函數在定義其實現方法時，可以采用可變參數列表varargin（即 variable-arguments-input），這可以使得無需指明可能的參數個數（更別提類型），在函數結構內部再進行識別和操作。當然，可以想見此時的varargin肯定不是普通的只能組織浮點數數據的向量；它必須能夠組織多種多樣的數據——元胞數組(cell)。

　　存放類的定義的文件一般和類同名(stressTensor.m)；在該文件內部需要指明類的屬性（成員變量），在“屬性(properties)”關鍵詞下完成（見上代碼）；指明類的方法（成員函數），在“方法(methods)”關鍵詞下完成（同見上代碼）。比如，求應力矢量、求正應力、求剪應力：

function stress = getStress( obj, direction )
    % 求應力矢量
    % 輸入一個參數：待求面法向量direction；
    % 返回：該面上應力矢量stress（行向量）
    direction = direction / norm(direction);
    if size(direction, 1) > 1
        stress = (obj.tensorMat * direction).';
    else
        stress = direction * obj.tensorMat;
    end
end

function stressNorm = getNormal( obj, vec )
    % 求正應力
    % 輸入一個參數：待求面法向量vec；
    % 返回：該面上正應力stressNorm（行向量）
    normVal = obj.getNormalVal(vec);
    stressNorm = normVal*vec/norm(vec);
end

function stressShear = getShear( obj, vec )
    % 求剪應力
    % 輸入一個參數：待求面法向量vec；
    % 返回：該面上剪應力矢量stressShear（行向量）
    vunit = vec/norm(vec);
    stress = obj.getStress(vec);
    stressVal = dot(stress, vunit);
    stressNorm = stressVal*vunit;
    stressShear = stress - stressNorm;
end

　　以及求主應力的方法：

function pAxis = principalAxis(obj)
    % 求主應力及其對應方向
    % 空輸入；返回一個3*2元胞數組，第i行為：第i個主應力-第i個主應力對應方向
    pAxis = cell(3, 2);
    A = obj.tensorMat;
    [eigvec, eigval] = eig(A);
    pStress = diag(eigval);
    for i = 1:3
        pAxis{i, 1} = pStress(i);
        pAxis{i, 2} = eigvec(:, i).';
    end
end

　　我們注意到在類的函數中有參數obj，這個詞表明該函數中需要調用類自身的其他函數或成員變量，這和python中但凡是從屬於類的方法（成員函數），在定義類的過程中必須要包含self這個參數所起的作用是相同的。C++中用法類似的是類自身指針this，但是在參數中並不需要包含。包含了obj參數以后就可以在function內部調用自身的屬性和方法了，所用的語法分別是obj.property和obj.method(arguments)。就我所知，MATLAB中類函數的定義中obj不是強求包含的，但是如果不包含，這個函數放在類中的意義就不明確了，除類內重載的函數外，往往可以放在類外（對比Python的靜態(static)函數）。

 

三、類函數的重載和運算符重載

　　最后，談到MATLAB就不能不談到數值計算，即使是面向對象編程也不例外；談到數值計算，就必須談到運算符重載，因為只有這樣才能真正地把自己寫的類融入到MATLAB的一套運算語法中去。MATLAB里的運算符重載比其他許多語言都顯得直截了當得多，因為它的每個運算符都對應一個函數，這個函數和普通的函數沒有什么區別。比如+運算符，對應的函數文件就是plus.m，函數名即plus。因此在自定義類中，運算符重載和函數重載別無二致，只要重寫一個位於這個類下的同名方法即可，例如：

function result = plus( obj1, obj2 )
    A = obj1.tensorMat + obj2.tensorMat;
    result = stressTensor(A);
end

　　只要確定plus是從屬於stressTensor的方法，就實現了+運算符重載。

　　具體有哪些運算符可以重載，對應的函數名是什么可以參見運算符重載參考頁。或者下面這個很大的表也提供了相應的信息：

 

四、檢驗和實際調用

　　接下來，只需要把以上這一系列函數放在類定義體(classdef stressTensor)的方法(methods)關鍵詞下就可以了。我們可以寫一個腳本調用：

import stressTensor.*;    % 首先，import這個類及其內容，可使用通配符.*導入類下的所有內容
T1 = stressTensor();    % 無參構造函數，參見以上構造函數定義時nargin==0的情形
A = 3*eye(3) + diag([1, 2], 1) + diag([1, 2], -1);
T2 = stressTensor(A);    % 含參構造函數
T = T1 + T2;    % 調用重載的加法
disp(T.tensorMat);    % 顯示此時的應力張量
%    輸出：
%     3     1     0
%     1     3     2
%     0     2     3
v = [1, 1, 1];
stress = T.getStress(v);    % 調用求應力矢量的方法
disp(stress);
%    輸出：
%    2.3094    3.4641    2.8868
disp(T.getShearVal(v));    % 調用求剪應力的方法
%    輸出：
%   -0.5774    0.5774   -0.0000　　

　　得到的結果符合預期。

 

五、類函數拆分為文件和組織方式

　　最后，關於文件的拆分。

　　一個類也許會很大，包含一些屬性和數量龐大和規模龐大的方法，這肯定會造成我們剛剛寫的stressTensor.m這樣儲存類的信息的文件越來越大、越來越亂，不便於管理。對此，C++的處理方法是，將類的定義和類函數的聲明放在頭文件里，函數的實現放在源文件里，而且可以通過多個源文件完成函數實現。MATLAB更簡單，所有的函數都可以原封不動地拆出來，把所有的代碼剪切到一個新文件里就行了，新文件的文件名設置為函數的名稱。但是！為了讓MATLAB編輯器明白它們都是從屬於一個類的函數，所有關於一個類的方法和這個類的定義本身需要放在一個文件夾里，這個文件夾的名字為“@"+類名，@是為了讓編譯器明白這是一個類的文件夾；比如之前的應力張量，文件夾名為"@stressTensor"。下圖就是我的這個stressTensor類的組織方式。