MATLAB算符重載和量綱分析

All comments and opinions expressed on Zhihu are mine alone and do not necessarily reflect those of my employers, past or present.
本文內容所有內容僅代表本人觀點，和Mathworks無關
工程科學計算中，一般物理量都是有單位的，比如：速度單位是米每秒LT^{-1}， 由長度基本量綱和時間基本量綱構成。加速度單位是米每秒平方，其量綱是LT^{-2}。 加速度乘以質量得到力，單位牛頓，其量綱是MLT^{-1}： F=ma.單位和量綱的運算遵從量綱法則：只有量綱相同的物理量，才能彼此相加、相減。也就是說我們不可以把速度和加速度相加：

?= v + an

工程科學計算中，如果不小心對不同單位的物理量做了加減運算，不但結果是錯誤的，而且應用到實際也可能會帶來危險。所以有必要建立這樣的一個單位和量綱系統，在計算的過程中可以攜帶單位，計算得到的結果也是有單位的物理量，並且在不小心對不同單位的物理做了加減運算時，該系統能夠終止計算並且提示錯誤。 具體的，這個量綱系統需要幫助我們實現以下的功能（偽碼）：

% 基本的需求偽代碼na0 = 2 ; % 變數a0表示值為2的加速度， 量綱是[LT^{-2}]nm1 = 3 ; % 變數m1表示值為1的質量， 量綱是[M]nm2 = 4 ; % 變數m2表示值為2的質量， 量綱是[M]nm3 = m1 + m2 ; % 允許相同量綱的物理量相加，結果的量綱仍然是[M]nf4 = a0 * m1 ; % 結果的單位是牛頓，量綱是[MLT^{-2}]nnx = a0 + m1 ; % 禁止，報錯ny = m1 - f4 ; % 禁止，報錯n

如何表示量綱

國際標準量綱制規定了物理量的基本量綱是：質量，長度，時間，電荷，溫度，密度和物質的量，其他物理量的量綱在這些基本量綱的基礎上複合而成。最簡單的量綱系統的可以用一個1 X 7的整數數組來表示：

% 初步設計nmass = [1,0,0,0,0,0,0];nlength = [0,1,0,0,0,0,0];ntime = [0,0,1,0,0,0,0];ncharge = [0,0,0,1,0,0,0];ntemperature = [0,0,0,0,1,0,0];nintensithy = [0,0,0,0,0,1,0];namount_of_substance = [0,0,0,0,0,0,1];n

而複合量綱可以用以上基本量綱的冪次來表示，比如：

% 複合量綱是基本量綱的冪積nveloctiy = [0,1,-1,0,0,0,0];nacceleration = [0,1,-2,0,0,0,0];nforce = [1,1,-2,0,0,0,0];n...n

如果我們需要讓每個物理量在計算中攜帶這些單位，最容易先想到就是用結構體表示這樣的物理量：

% 基本的需求偽代碼na0.value = 2 ; % 值na0.dim = [0,1,-2,0,0,0,0]; % 單位是加速度nnm1 = 3 ; % 值nm1.dim = [1,0,0,0,0,0,0]; % 單位是是質量n

但是MATLAB中，結構體不可以直接用來做代數運算，如果嘗試如下的運算，MATLAB將報錯:

% 結構體沒有重載算符n>> f4 = m1 * a0n Undefined operator * for input arguments of type struct.n

於是不可避免的，我們還需要設計如下的運算函數來處理結構體之間的算術運算

% 結構體無法重載算符，所以需要用函數來實現相互的算術運算nnew_s = dim_plus(s1,s2); % 結構體之間的加nnew_s = dim_minus(s1,s2); % 結構體之間的減nnew_s = dim_product(s1,s2;) % 結構體之間的乘nnew_s = dim_divide(s1,s2); % 結構體之間的除n

大致的，dim_plus函數可以這樣實現：函數內部首先檢查作為輸入的兩個的結構體的dim域，比較它們是否相同，以此達到強制量綱運算的法則的目的：

% dim_plusnfunction new_s = dim_plus(s1,s2)n % .... 省略其它對輸入的檢查n assert(s1.dim == s2.dim) % 強制單位之間的加法只能在相同量綱之間進行n new_s.dim = s1.dim; % 結果的量綱不變n new_s.value = s1.value + s2.value;nendn

使用如下：

% 如果在不同的單位之間做加法運算 n>> s = dim_plus(a0,m1)n Error using assertn The condition input argument must be a scalar logical.n Error in dim_plus (line 2)n assert(s1.dim == s2.dim)n

結構體來表示量綱的缺點是顯而易見的：直接的算術運算被替換成了第1行的dim_plus的函數調用，閱讀和使用都很麻煩。回到面向對象的程序設計中來，我們已經知道：封裝對數據的運算恰好是面向對象程序設計的強項，而且面向對象還提供了對算符的重載功能，使得我們可以自己設計如何解釋形如 m1 + m2 , f = m1 * a0 這樣的運算，所以我們其實可以利用面向對象來設計量綱系統。

需求和設計:加法和減法

在本章的後續小節中，我們將採用這樣的一個工作流程：首先描述程序的需求（requirement），即從外部看上去這個量綱系統的工作方式(也叫做Functional API設計) ，再利用需求去引導設計。首先我們需要一個量綱類，該類應該接受 1 X 7的數組作為輸入，表示該對象的實際量綱：

% 設計Construct允許的輸入n>> Dimension([1,0,0,0,0,0,0]); % massn>> Dimension([0,1,0,0,0,0,0]); % lengthn>> Dimension([0,0,1,0,0,0,0]); % timen

並且我們規定，該只接受1X7的行向量，不接受元胞，不接受列向量：

%設計Constructor禁止的輸入n>> Dimension({1,2}) % 不接受元胞n>> Dimension([2,2,2,2,2,2,2,2,2]) % 不接受其它大小的行向量n>> Dimension([0;0;1;0;0;0;0]); % 不接受列向量n

滿足上麵條件的Dimension類如下，其中第7行使用了validateattributes函數來驗證函數的輸入必須是1 X 7的數組，

classdef Dimensionn propertiesn value % dimension valuen endn methodsn function obj = Dimension(input)n validateattributes(input,{numeric},{size,[1 7]});n obj.value = input;n endn endnendn

Dimension類的重要功能就是對算術運算的進行單位檢查：比如加減運算只能在相同量綱的物理量之間進行，要滿足的所有功能的清單如下：

%加減運算nt1 = Dimension([0,0,1,0,0,0,0]) ; % timent2 = Dimension([0,0,1,0,0,0,0]) ; % timent3 = t1 + t2; % 允許，結果量綱是 [0,0,1,0,0,0,0]nnm1 = Dimension([1,0,0,0,0,0,0]) ; % massnm2 = Dimension([1,0,0,0,0,0,0]) ; % massnm3 = m1 + m2; % 允許，結果量綱是 [1,0,0,0,0,0,0]nnnt1 = Dimension([0,0,0,0,0,0,0]) ; % scalarnl1 = Dimension([0,1,0,0,0,0,0]) ; % lengthnx = t1 + l1; % 禁止，報錯nnl1 = Dimension([0,1,0,0,0,0,0]) ; % lengthnm1 = Dimension([1,0,0,0,0,0,0]) ; % massnx = l1 - m2; % 禁止，報錯n

為了達到上述的需求，我們必須重載Dimension類的plus和minus算符：

classdef Dimensionn% value object, can do compare directlyn propertiesn value % dimension valuen endn methodsn .... 其餘略n function newunit = plus(o1,o2)n isequalassert(o1,o2);n newunit = o1; % 結果的量綱等於o1的量綱或者o2的量綱n endnn function newunit = minus(o1,o2)n isequalassert(o1,o2);n newunit = o1; % 結果的量綱等於o1的量綱或者o2的量綱n endnendnn% 類Dimension的局部函數nfunction isequalassert(v1,v2)n if isequal(v1,v2)nelsen error(Dimension:DimensionMustBeTheSame,);nendnnendn

說明如下

• 第1行中，我們把該Dimension類設計成了Value類，這是為了方便第21行直接對兩個輸入的量綱進行比較.(handle類對象之間的比較具有不同的意義)
• 第9和第14行中，我們調用了類的局部的函數(請參考相關章節) isequalassert，在每次進行計算之前，檢查運算的量綱是否相同，如果不同，則報錯。
• 量綱的加減運算，結果量綱不變，所以第10，15行直接構造一個新的Dimension對象，等於原來的量綱值，作為結果返回。

目前為止，Dimension類只是一個表示單位的類，不包括物理量的實際值。上面的設計加法和減法也僅僅限於單位之間的加減法。為了表示物理量的實際的值，我們還需要一個Quantity類，該類既包括物理量的值，也包括物理量的Dimension。該類需要滿足如下簡單構造：提供value和量綱對象做為構造函數的輸入，返回Quantity對象

m1 = Quantity(2,Dimension([1,0,0,0,0,0,0])); % 變數q1表示2千克的質量na1 = Quantity(3,Dimension([0,1,-2,0,0,0,0]));% 變數a1表示加速度為3n

我們應該可以對Quantity對象進行加減法，如下所示：

m1 = Quantity(2,Dimension([1,0,0,0,0,0,0])); % 變數m1表示2千克的質量nm2 = Quantity(4,Dimension([1,0,0,0,0,0,0])); % 變數m2表示4千克的質量nm3 = m1 + m2 ; % 結果m3表示6千克的質量n

在Dimension類的設計的基礎上，容易寫出Quantity類的：構造函數和加法減法的重載：

classdef Quantityn propertiesn valuen unitn endn methodsn function obj = Quantity(value,unit)n obj.value = value;n obj.unit = unit;n endn function results = plus(o1,o2)n results = Quantity(o1.value+o2.value,o1.unit+o2.unit);n endn function results = minus(o1,o2)n results = Quantity(o1.value-o2.value,o1.unit-o2.unit);n endn endnendn

其中第12行和第15行的第2個參數把對Quantity的單位的加減法的計算轉到了Dimension的加減法計算函數上去。

需求和設計:乘法和除法

量綱運演算法則對乘除法的規定是：複合量綱是其它的量綱的冪積，這也就是說：

% Dimension對象的乘除法需求ns1 = Dimension([0,1,0,0,0,0,0]) ; % lengthnt1 = Dimension([0,0,1,0,0,0,0]) ; % timenv1 = s1/t1 ; % 結果是速度 量綱是[0,1,-1,0,0,0,0]nns2 = v1 * t1 ; % 結果是長度 量綱是[1,0,0,0,0,0,0]nnm1 = Dimension([1,0,0,0,0,0,0]) ; % massna1 = Dimension([0,1,-2,0,0,0,0]) ; % accelerationnf1 = m1 * a1 ; % 結果單位是力 量綱是[1,1,-2,0,0,0,0]n

上述是簡單的量綱本身的乘除法，而實際的物理量Quantity對象的乘除法，需要滿足如下需求：

% Qauntity對象的乘除法的需求ns1 = Quantity(2, Dimension([0,1,0,0,0,0,0])) ; % 2米nt1 = Quantity(4, Dimension([0,0,1,0,0,0,0])) ; % 4秒nv1 = s1/t1 ; % 結果Qunatity對象 值0.5 單位是速度 量綱[0,1,-1,0,0,0,0]nns2 = v1 * t1 ; % 結果Qunatity對象 值2 單位是長度 量綱[1,0,0,0,0,0,0]nnm1 = Quantity(3,Dimension([1,0,0,0,0,0,0])) ; % 3千克na1 = Quantity(5,Dimension([0,1,-2,0,0,0,0])) ; % 5米每秒平方nf1 = m1 * a1 ; % 結果Qunatity對象 值15 單位是牛 量綱[1,1,-2,0,0,0,0]n

除了可以用Quantity的構造函數來生產Quantity對象，根據我們計算的書寫習慣，我們還期望Quantity類還支持其它的對象的方法。比如8秒鐘這個物理量，可以用Quantity類的構造函數來構造：

t1 = Quantity(8, Dimension([0,0,1,0,0,0,0]))n

當然我們還知道8秒鐘，其實還可以通過 4秒鐘 $times 2$ ,或者 8 $times$ 一個單位時間秒來表示，這裡的2和8都是標量，或者普通的數字。該需求反映在代碼上，就是要求下面的計算都能返回value為8，unit是Dimension([0,0,1,0,0,0,0]) 的Quantity對象

% Quantity類對象和scalar的乘法得到Quantity對象nt1 = Quantity(2,Dimension([0,0,1,0,0,0,0]) * 4nt2 = 4 * Quantity(2,Dimension([0,0,1,0,0,0,0])nn% Dimension類對象和scalar的乘法得到Quantity對象nt3 = 8 * Dimension([0,0,1,0,0,0,0]nt4 = Dimension([0,0,1,0,0,0,0] * 8n

目前為止我們對這個量綱系統的需求描述完畢，根據這些需求，我們給Dimension的類新添了第21到48行，其 完整設計如下：

classdef Dimensionn % value object, can do compare directlyn propertiesn value % dimension valuen endn methodsn function obj = Dimension(input)n validateattributes(input,{numeric},{size,[1 7]});n obj.value = input;n endnn function newunit = plus(o1,o2)n isequalassert(o1,o2);n newunit = o1;n endnn function newunit = minus(o1,o2)n isequalassert(o1,o2);n newunit = o1;n endnn function newObj = mtimes(o1,o2)n validateattributes(o1,{numeric,Dimension},{});n validateattributes(o2,{numeric,Dimension},{});n if isnumeric(o1) && isa(o2,Dimension)n newObj = Quantity(o1,o2);n elseif isnumeric(o2) && isa(o1,Dimension)n newObj = Quantity(o2,o1) ;n elseif isa(o1,Dimension) && ...n isa(o2,Dimension)n newObj = Dimension(o1.value + o2.value);n elsen % will not reach here, due to dispatch rulesn endn endnn function newObj = mrdivide(o1,o2)n validateattributes(o1,{numeric,Dimension},{});n validateattributes(o2,{numeric,Dimension},{});n if isnumeric(o1) && isa(o2,Dimension)n newObj = Quantity(o1,Dimension(-o2.value));n elseif isnumeric(o2) && isa(o1,Dimension)n newObj = Quantity(1/o2,o1) ;n elseif isa(o1,Dimension) ...n && isa(o2,Dimension)n newObj = Dimension(o1.value - o2.value);n elsen % will not reach here, due to dispatch rulesn endn endn endnendnnfunction isequalassert(v1,v2)n if isequal(v1,v2)n elsen error(Dimension:DimensionMustBeTheSame,);n endnendn

Dimension類設計說明如下

• 第22，23，36，37行確保做乘除運算時，運算符要麼是兩個Dimension對象

v1 = s1/t1 ;ns2 = v1 * t1 ;n

要麼是Dimension對象和簡單的標量

t3 = 8 * Dimension([0,0,1,0,0,0,0]nt4 = Dimension([0,0,1,0,0,0,0] * 8n

根據上述需求，我們給Quantity的類新添了17-45行，其完整設計如下：

classdef Quantityn propertiesn valuen unitn endn methodsn function obj = Quantity(value,unit)n obj.value = value;n obj.unit = unit;n endn function results = plus(o1,o2)n results = Quantity(o1.value+o2.value,o1.unit+o2.unit);n endn function results = minus(o1,o2)n results = Quantity(o1.value-o2.value,o1.unit-o2.unit);n endn function newObj = mtimes(o1,o2)n [o1,o2] = converter_helper(o1,o2);n newObj = Quantity(o1.value*o2.value,o1.unit*o2.unit);n endn function newObj = mrdivide(o1,o2)n [o1,o2] = converter_helper(o1,o2);n newObj = Quantity(o1.value/o2.value,o1.unit/o2.unit);n endn endnendnn% convert input into Quantity objectnfunction [o1,o2] = converter_helper(o1,o2)n validateattributes(o1,{numeric,Quantity,Dimension},{});n validateattributes(o2,{numeric,Quantity,Dimension},{});nn% convert numeric input into Quantity objectn if isnumeric(o1)n o1 = Quantity(o1,Dimension([0,0,0,0,0,0,0]));n elseif isnumeric(o2)n o2= Quantity(o2,Dimension([0,0,0,0,0,0,0]));n endnn% convert dimension input into Quantity objectn if isa(o1,Dimension)n o1 = Quantity(1,o1);n elseif isa(o2,Dimension)n o2= Quantity(1,o2);n endnendn

Quantity類說明如下

• 第18，22行對輸入進行預處理，調用局部函數converter_helper確保兩個參數們統統都是Quanit對象，再分別對其value和unit部分做計算。第19, 23行返回新的Quantity對象。
• 局部函數converter_helper從第29行開始，第30，31行validateattributes限制操作數只能是簡單的數字,Quantity或者Dimension對象。如果有任何一個輸入是簡單的數字, 第34到38行把它轉成Quantity對象，量綱為Scalar。如果有任何一個輸入是Dimension對象，第41到45行把它轉成Quantity對象，其值是1，其Dimension不變.

在單元測試和以測試為驅動的章節中，我們還將重新使用這個量綱系統，在那裡，我們將重新審視本節的工作流程，從而總結一套更加規範的工業設計和開發流程。
推薦閱讀：