MATLAB 超程式設計介紹

這篇文章對 Matlab 中的超程式設計進行了簡單的介紹。Matlab 是一個古老而又高度專業化的語言。由於這一原因，缺乏很多在現代或者通用語言中擁有的特性。然而，用一些簡單的工具，我們可以發現 Matlab 也可以足夠靈活去進行非常簡單的超程式設計。

什麼是超程式設計？為什麼用 Matlab 來做

粗淺的說，超程式設計是將程式視為資料的過程——意味著一個程式可以像一個普通的資料片段一樣被製造並且被操作。像 Ruby 語言就顯示出了一些作為超程式設計語言的特性。Matlab 未能提供這樣的一套工具用於進行超程式設計。然而事實上用一點創造性就能實現，用一些簡單的工具箱完成這一點。

注意如果把這一部分的內容稱為超程式設計的話，可能會有一些誇張，但是我認為內容上已經足夠接近去自證這樣的標籤。

前提條件

本文的程式碼能執行是基於如下條件的:

• 函式是可例項化的（first class objects）
• 函式可以被匿名定義
• 符號表示式可以變作程式執行.

粗略的講， 這些使我們可以定義符號表示式，修改他們，而且使他們變為函式。這就是這篇文章的核心觀點

使用和修改符號表示式

首先定義一個符號變數

x = sym('x');

MATLAB 將會將 x 當做符號來處理，而不是數字。 然後我們可以構造更復雜的表示式，例如:

>> x + x
ans =
2*x
>> x * x
ans =
x^2
>> sin(x)^2 + cos(x)^2
ans =
sin(x)^2 + cos(x)^2

當然，最後一個表示式恆等於1. MATLAB 也知道:

>> simplify(sin(x)^2 + cos(x)^2)
ans =
1

而且可以做微積分:

>> diff(3*x^2 + sin(x))
ans =
6*x + cos(x)
>> int(6*x + cos(x))
ans =
sin(x) + 3*x^2
>> taylor(exp(x), x, 0, 'Order', 10)
ans =
x^9/362880 + x^8/40320 + x^7/5040 + x^6/720 + x^5/120 + x^4/24 + x^3/6 + x^2/2 + x + 1

symbolic 工具箱有一個叫 matlabFunction 的函式，可以把表示式專為一個函式:

>> p = matlabFunction(x^2 + 2*x + sin(x))
p =
@(x)x.*2.0+sin(x)+x.^2
>> p(0.2)
ans =
0.6387

有此利器在手，我們就可以編寫，能寫程式的程式了。

傅立葉級數逼近

傅立葉級數是一種函式逼近，由不同週期的正弦、餘弦進行線性組合構成。 一個實際函式 f(x) 可以由如下級數逼近:

f(x)≈a02+∑n=1∞ancos(nx)+∑n=1∞bnsin(nx)

where:

a0=1π∫π−πf(x)dx an=1π∫π−πf(x)cos(nx)dx bn=1π∫π−πf(x)sin(nx)dx

這個近似已經精確地由這幾個函式描述了， 但是我們準備寫一個程式來近似一個隨機的真實的函式，就用這個級數的截斷的版本。

超程式設計實現傅立葉級數近似

現在，我們要實現一個函式，接收一個函式作為引數，返回另一個函式，是輸入函式的截斷的傅立葉級數近似。 有些更好的方法來實現超程式設計， 這裡只是作為超程式設計的優雅的示例。

這個函式應按照如下形式工作

>> g = fourierapprox(f, N);

其中 f 是一個函式， N 是 f逼近中的最高項，g 為得到的近似函式。

首先， 我們建立一個符號表示式（最終會被轉為MATLAB函式）， 定義為一個獨立變數x 。然後新增series的第一項:

function [fn] = fourierapprox(f, N)
    series = sym();
    x = sym('x');
    a_0 = (1/pi) * integral(f, -pi, pi);
    series = series + (1/2)*a_0;
    ...
end

(注意，我們使用 MATLAB自帶的數值方法來計算必要的積分。)

下一步， 我們需要向series中新增有限項 (N)的正弦跟餘弦的線性組合 :

function [fn] = fourierapprox(f, N)
    series = sym();
    x = sym('x');
    a_0 = (1/pi) * integral(f, -pi, pi);
    series = series + (1/2)*a_0;
    for n = 1:N
        a_n = (1/pi) * integral(@(x) f(x) .* cos(n*x), -pi, pi);
        series = series + a_n * cos(n*x);
        b_n = (1/pi) * integral(@(x) f(x) .* sin(n*x), -pi, pi);
        series = series + b_n * sin(n*x);
    end
    fn = matlabFunction(series);
end

最後一行很關鍵: 他將一個單純的符號變數轉為一個可執行的函式。 這意味著 fn 表現的跟硬編碼的函式一樣! 稱之為 ‘超程式設計’ 是沒問題的， 因為我們通過程式生成了這個函式。

執行這個函式

我們來用剛才的函式生成一些近似函式看看。

下面這個函式將使用特定的間距在同一座標系內繪出原函式跟傅立葉近似函式

function [ output_args ] = plotfourierapprox(f, N, a, b)
    fplot(f, [a b], 'r');
    hold on;
    pn = fourierapprox(f, N);
    fplot(pn, [a b], 'b');
end

方波

>> plotfourierapprox(@square, 5, -5, 5)

>> plotfourierapprox(@square, 15, -5, 5)

f(x)=x3

>> plotfourierapprox(@(x) x.^3, 10, -4, 4)

最後的話

超程式設計不是MATLAB中典型的程式設計正規化。 我想展示的是，這不但可能，而且實用。