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

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

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

　　在彈性力學裏，應力張量是描述彈性體內部某一點的應力情況的；它和過該點且法向量爲 $\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()這樣的方法調用豈不是美滋滋？這就是要學習——

 

 

2、定義應力張量的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)函數）。

 

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

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

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

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

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

 

4、檢驗和實際調用

　　接下來，只須要把以上這一系列函數放在類定義體(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　　

　　獲得的結果符合預期。

 

5、類函數拆分爲文件和組織方式

　　最後，關於文件的拆分。

　　一個類也許會很大，包含一些屬性和數量龐大和規模龐大的方法，這確定會形成咱們剛剛寫的stressTensor.m這樣儲存類的信息的文件愈來愈大、愈來愈亂，不便於管理。對此，C++的處理方法是，將類的定義和類函數的聲明放在頭文件裏，函數的實現放在源文件裏，並且能夠經過多個源文件完成函數實現。MATLAB更簡單，全部的函數均可以原封不動地拆出來，把全部的代碼剪切到一個新文件裏就好了，新文件的文件名設置爲函數的名稱。可是！爲了讓MATLAB編輯器明白它們都是從屬於一個類的函數，全部關於一個類的方法和這個類的定義自己須要放在一個文件夾裏，這個文件夾的名字爲「@"+類名，@是爲了讓編譯器明白這是一個類的文件夾；好比以前的應力張量，文件夾名爲"@stressTensor"。下圖就是個人這個stressTensor類的組織方式。