MATLAB 元编程介绍

破旧瓷碗

这篇文章对 Matlab 中的元编程进行了简单的介绍。Matlab 是一个古老而又高度专业化的语言。由于这一原因，缺乏很多在现代或者通用语言中拥有的特性。然而，用一些简单的工具，我们可以发现 Matlab 也可以足够灵活去进行非常简单的元编程。

什么是元编程？为什么用 Matlab 来做

粗浅的说，元编程是将程序视为数据的过程——意味着一个程序可以像一个普通的数据片段一样被制造并且被操作。像 Ruby 语言就显示出了一些作为元编程语言的特性。Matlab 未能提供这样的一套工具用于进行元编程。然而事实上用一点创造性就能实现，用一些简单的工具箱完成这一点。

注意如果把这一部分的内容称为元编程的话，可能会有一些夸张，但是我认为内容上已经足够接近去自证这样的标签。

panlingrui
翻译于 2015/05/18 09:23

其它翻译版本 (1)
前提条件
本文的代码能运行是基于如下条件的:

函数是可实例化的（first class objects）

函数可以被匿名定义

符号表达式可以变作程序执行.

粗略的讲， 这些使我们可以定义符号表达式，修改他们，而且使他们变为函数。这就是这篇文章的核心观点

使用和修改符号表达式
首先定义一个符号变量

1. x = sym('x');

复制代码

MATLAB 将会将 x 当做符号来处理，而不是数字。 然后我们可以构造更复杂的表达式，例如:

1. >> x + x
   
2. ans =
   
3. 2*x
   
4. >> x * x
   
5. ans =
   
6. x^2
   
7. >> sin(x)^2 + cos(x)^2
   
8. ans =
   
9. sin(x)^2 + cos(x)^2

复制代码

当然，最后一个表达式恒等于1. MATLAB 也知道:

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

复制代码

而且可以做微积分:

1. >> diff(3*x^2 + sin(x))
   
2. ans =
   
3. 6*x + cos(x)
   
4. >> int(6*x + cos(x))
   
5. ans =
   
6. sin(x) + 3*x^2
   
7. >> taylor(exp(x), x, 0, 'Order', 10)
   
8. ans =
   
9. 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 的函数，可以把表达式专为一个函数:

1. >> p = matlabFunction(x^2 + 2*x + sin(x))
   
2. p =
   
3. @(x)x.*2.0+sin(x)+x.^2
   
4. >> p(0.2)
   
5. ans =
   
6. 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

这个近似已经精确地由这几个函数描述了， 但是我们准备写一个程序来近似一个随机的真实的函数，就用这个级数的截断的版本。

元编程实现傅里叶级数近似
现在，我们要实现一个函数，接收一个函数作为参数，返回另一个函数，是输入函数的截断的傅里叶级数近似。 有些更好的方法来实现元编程， 这里只是作为元编程的优雅的示例。

这个函数应按照如下形式工作

1. >> g = fourierapprox(f, N);

复制代码

其中 f 是一个函数， N 是 f逼近中的最高项，g 为得到的近似函数。
首先， 我们创建一个符号表达式（最终会被转为MATLAB函数）， 定义为一个独立变量x 。然后添加series的第一项:

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

复制代码

(注意，我们使用 MATLAB自带的数值方法来计算必要的积分。)

下一步， 我们需要向series中添加有限项 (N)的正弦跟余弦的线性组合 :

1. function [fn] = fourierapprox(f, N)
   
2.     series = sym();
   
3.     x = sym('x');
   
4.     a_0 = (1/pi) * integral(f, -pi, pi);
   
5.     series = series + (1/2)*a_0;
   
6.     for n = 1:N
   
7.         a_n = (1/pi) * integral(@(x) f(x) .* cos(n*x), -pi, pi);
   
8.         series = series + a_n * cos(n*x);
   
9.         b_n = (1/pi) * integral(@(x) f(x) .* sin(n*x), -pi, pi);
   
10.         series = series + b_n * sin(n*x);
    
11.     end
    
12.     fn = matlabFunction(series);
    
13. end

复制代码

最后一行很关键: 他将一个单纯的符号变量转为一个可执行的函数。 这意味着 fn 表现的跟硬编码的函数一样! 称之为 ‘元编程’ 是没问题的， 因为我们通过程序生成了这个函数。

执行这个函数
我们来用刚才的函数生成一些近似函数看看。

下面这个函数将使用特定的间距在同一坐标系内绘出原函数跟傅里叶近似函数

1. function [ output_args ] = plotfourierapprox(f, N, a, b)
   
2.     fplot(f, [a b], 'r');
   
3.     hold on;
   
4.     pn = fourierapprox(f, N);
   
5.     fplot(pn, [a b], 'b');
   
6. end

复制代码

方波

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

复制代码

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

复制代码

f(x)=x3

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

复制代码