Octave笔记

Octave可以看作是开放源代码版本的Matlab，它的命令几乎和Matlab完全一样，因此可以用Matlab的书籍和文档来学习Octave。另外，绝大多数Matlab的代码也可以不经修改在Octave中运行。对于初学者来说，下面这些笔记可以帮助快速上手Octave。更多的功能，可以参考Octave文档或者其他Matlab的资料。

1. 脚本文件，以下面一行开始：

   #!/usr/bin/octave -qf

2. 注释，以%或者#开始
3. 虚数单位：i, j, I, J都可以 两种方法定义复数：3 + 4i或者complex(3, 4)
4. 定义矩阵：

   octave:1> a = [1, 2; 3, 4]
   a =
   
   1   2
   3   4

5. 随机数：rand()或者rand(m, n)
6. 定义一段连续序列：

   octave:1> 1:5
   ans =
   
   1   2   3   4   5

   还可以指定步长：

   octave:2> 1:2:10
   ans =
   
   1    3    5    7    9

7. 字符串：单引号或者双引号
8. 定义Data Structure：

   octave:1> x.name = ”Yusuf”
   octave:2> x.age = 33
   octave:3> x.phone = “46103”
   octave:4> x
   x =
   {
   name = Yusuf
   age = 33
   phone = 46103
   }

9. 定义Cell Array：

   octave1:> a = {”random matrix”, rand(2)}
   a =
   
   {
   [1,1] = random matrix
   [1,2] =
   
   0.417140   0.460380
   0.337105   0.078835
   
   }
   
   octave2:> a{2}(1,2)
   ans =  0.460380

10. 声明全局变量：

    global c = 3e8

11. 算术符号

    x + y	加法
    x - y	减法
    x * y	矩阵乘法
    x .* y	元素相乘
    x / y	右除，相当于(inverse (y’) * x’)’
    x ./ y	元素右除
    x \ y	左除，相当于inverse (x) * y
    x .\ y	元素左除，每个y元素被相应的x元素除
    x ^ y 或 x ** y	乘方
    x .^ y 或 x .** y	元素乘方
    x’	转置复共轭
    x.’	转置

12. 比较符号

    x < y	小于
    x <= y	小于等于
    x == y	等于
    x >= y	大于等于
    x > y	大于
    x != y 或 x ~= y	不等于

13. if语句
• 形式一：

  if (condition)
    then-body
  endif

• 形式二：

  if (condition)
    then-body
  else
    else-body
  endif

• 形式三：

  if (condition)
    then-body
  elseif (condition)
    elseif-body
  else
    else-body
  endif

14. switch语句

    switch expression
      case label
        command_list
      case label
        command_list
      …
      otherwise
        command_list
    endswitch

15. while语句

    while (condition)
      body
    endwhile

16. do-until语句

    do
      body
    until (condition)

17. for语句

    for var = expression
      body
    endfor

18. 定义函数
• 形式一：无自变量 function name body endfunction
• 形式二：带自变量 function name (arg-list) body endfunction
• 形式三：单个返回值 function ret-var = name (arg-list) body endfunction
• 形式四：多个返回值 function [ret-list] = name (arg-list) body endfunction
19. 可变自变量长度函数 function val = smallest (varargin) val = min ([varargin{:}]); endfunction
20. 函数Handles（指向函数的指针） @function-name
21. 匿名函数 @(argument-list) expression
22. 终端输出 disp (x)
23. 写读包含分隔符的数据文件 dlmwrite (file, a) dlmwrite (file, a, delim, r, c) dlmwrite (file, a, key, val ...) dlmwrite (file, a, "-append", ...) a是待写矩阵，delim是分隔符，默认为逗号。r, c分别为起始的空行和空列。key包含下列参数：
• "append" 即上面的"-append"，"on"或者"off"
• "delimiter" 即上面的delim
• "newline" 换行符
• "roffset" 即上面的r
• "coffset" 即上面的c
• "precision" 精确度，可以用类似fprintf的格式字符串，或者有效数字位数
data = dlmread (file) data = dlmread (file, sep, range) range为四元素的矢量：[R0, C0, R1, C1]，为数据块左上角和右下角的标记，从0开始
24. Octave可以使用大多数类似C的I/O函数（fprintf，fopen等）
25. 常用数学函数

    exp (x)	指数函数	sinh (x)	双曲正弦
    log (x)	自然对数	cosh (x)	双曲余弦
    log10 (x)	以10为底的对数	tanh (x)	双曲正切
    log2 (x)	以2为底的对数	coth (x)	双曲余切
    sqrt (x)	平方根	sech (x)	双曲正割
    abs (z)	绝对值	csch (x)	双曲余割
    arg (z)	辐角	asinh (x)	反双曲正弦
    conj (z)	复共轭	acosh (x)	反双曲余弦
    real (z)	复数实部	atanh (x)	反双曲正切
    imag (z)	复数虚部	acoth (x)	反双曲余切
    sin (x)	正弦	asech (x)	反双曲正割
    cos (x)	余弦	acsch (x)	反双曲余割
    tan (x)	正切	sum (x)	求和
    cot (x)	余切	prod (x)	求积u
    sec (x)	正割	ceil (x)	不小于x的最小整数
    csc (x)	余割	floor (x)	不大于x的最大整数
    asin (x)	反正弦	factorial (n)	阶乘
    acos (x)	反余弦	max (x)	最大值
    atan (x)	反正切	min (x)	最小值
    acot (x)	反余切	round (x)	四舍五入
    asec (x)	反正割	sign (x)	符号函数
    acsc (x)	反余割
    atan2 (y, x)	计算atan (y / x)
    sind (x),…	以度数为单位的三角函数

26. 特殊函数

    [a, ierr] = airy (k, z, opt)	Airy函数
    [j, ierr] = besselj (alpha, x, opt)	第一类Bessel函数
    [y, ierr] = bessely (alpha, x, opt)	第二类Bessel函数
    [i, ierr] = besseli (alpha, x, opt)	变形第一类Bessel函数
    [k, ierr] = besselk (alpha, x, opt)	变形第二类Bessel函数
    [h, ierr] = besselh (alpha, k, x, opt)	第一类（k=1）和第二类（k=2）Hankel函数
    beta (a, b)	Beta函数
    gamma (z)	Gamma函数
    bincoeff (n, k)	二项式系数
    erf (z)	误差函数
    legendre (n, x, normalization)	n阶Legendre函数，m = 0…n

27. 线性代数

    det (a)	行列式
    eig (a)	本征值
    dot (x, y)	点乘
    inv (a)	逆矩阵

28. 解非线性方程 [x, fval, info] = fsolve (fcn, x0, options) fcn是包含待解方程组的函数，x0是自变量初始值，x为返回的数值解，fval为返回的函数值，如果结果收敛，info为1，否则为其他值。例子： function y = f (x) y(1) = -2*x(1)^2 + 3*x(1)*x(2) + 4*sin(x(2)) - 6; y(2) = 3*x(1)^2 - 2*x(1)*x(2)^2 + 3*cos(x(1)) + 4; endfunction [x, fval, info] = fsolve (@f, [1; 2]) 结果： x = 0.57983 2.54621 fval = 6.1872e-08 -3.2708e-07 info = 1
29. 数值积分（各函数对应不同算法） [v, ier, nfun, err] = quad (f, a, b, tol, sing) quadl (f, a, b, tol) quadgk (f, a, b, prop, val, ...) quadl (f, a, b, tol)

统计笔记（一）

1. 协方差（covariance）定义：
   

   ${\displaystyle \mathrm{Cov}(Y_1,Y_2)=E[(Y_1-{\mu }_1)(Y_2-{\mu }_2)]}$

   
   其中 $E(\hspace{0.5em})$表示期望值， ${\mu }_1=E(Y_1)$， ${\mu }_2=E(Y_2)$。
2. 相关系数（coefficient of correlation）定义：
   

   ${\displaystyle \rho =\frac{\mathrm{Cov}(Y_1,Y_2)}{{\sigma }_1{\sigma }_2}}$

   
   其中 ${\sigma }_1$和 ${\sigma }_2$分别是 $Y_1$和 $Y_2$的标准差。
3. 任意分布样本的平均值
   $Y_1$， $Y_2$，...， $Y_n$是 $E(Y_i)=\mu$， $V(Y_i)={\sigma }^2$的独立随机变量，则平均值
   

   ${\displaystyle \overline{Y}=\frac{1}{n}\sum _{i=1}^n{Y}_i}$

   
   的期望值 $E(\overline{Y})=\mu$，方差 $V(\overline{Y})={\sigma }^2/n$。
4. 正态分布样本的平均值
   $Y_1$， $Y_2$，...， $Y_n$是来自正态分布（期望值为 $\mu$，方差为 ${\sigma }^2$）的样本，其平均值
   

   ${\displaystyle \overline{Y}=\frac{1}{n}\sum _{i=1}^n{Y}_i}$

   
   也呈正态分布，期望值 ${\mu }_{\overline{Y}}=\mu$，方差 ${\sigma }_{\overline{Y}}^2={\sigma }^2/n$。
5. $Y_1$， $Y_2$，...， $Y_n$同上，那么
   

   ${\displaystyle \sum _{i=1}^n{Z}_i^2=\sum _{i=1}^n{\left(\frac{Y_i-\mu }{\sigma }\right)}^2}$

   
   呈自由度为 $n$的 ${\chi }^2$分布。
6. $Y_1$， $Y_2$，...， $Y_n$同上，那么
   

   ${\displaystyle \frac{(n-1)S^2}{{\sigma }^2}=\frac{1}{{\sigma }^2}\sum _{i=1}^n(Y_i-\overline{Y}{)}^2}$

   
   呈自由度为 $n-1$的 ${\chi }^2$分布。其中，
   

   ${\displaystyle S^2=\frac{1}{n-1}\sum _{i=1}^n(Y_i-\overline{Y}{)}^2}$

   
   是样本方差。 $S^2$和 $\overline{Y}$为独立随机变量。
7. 学生 $t$分布(简称 $t$分布)定义：
   如果 $Z$是正态分布随机变量， $W$是自由度为 $\nu$的 ${\chi }^2$分布随机变量，如果 $Z$和 $W$互相独立，则
   

   ${\displaystyle T=\frac{Z}{\sqrt{W/\nu }}}$

   
   呈自由度为 $\nu$的 $t$分布。
8. 在上面的定义中，如果让 $Z=\sqrt{n}(\overline{Y}-\mu )/\sigma$， $W=(n-1)S^2/{\sigma }^2$，则由上面第4和第6条结论可以得出
   

   ${\displaystyle \sqrt{n}\left(\frac{\overline{Y}-\mu }{S}\right)}$

   
   呈自由度为 $n-1$的 $t$分布。
9. $F$分布定义：
   如果 $W_1$和 $W_2$为独立 ${\chi }^2$分布变量，自由度分别为 ${\nu }_1$和 ${\nu }_2$，则
   

   ${\displaystyle F=\frac{W_1/{\nu }_1}{W_2/{\nu }_2}}$

   
   被称作具有 ${\nu }_1$分子自由度和 ${\nu }_2$分母自由度的 $F$分布。
10. 中心极限定理（Central Limit Theorem）
    $Y_1$， $Y_2$，...， $Y_n$是独立的具有同样分布的随机变量，期望值 $E(Y_i)=\mu$，方差 $V(Y_i)={\sigma }^2<\infty$。定义
    

    ${\displaystyle U_n=\sqrt{n}\left(\frac{\overline{Y}-\mu }{\sigma }\right)}$

    
    其中
    

    ${\displaystyle \overline{Y}=\frac{1}{n}\sum _{i=1}^n{Y}_i}$

    
    那么，当 $n\to \infty$时， $U_n$趋近于正态分布。

让blogger支持MathML公式输出

我想在我的blog里显示公式怎么办？我可以把每个公式变成图片文件，然后把公式作为图片显示。我不需要截屏软件，OpenOffice可以很方便的把包含公式的文档转换为包含图片的HTML文件。但是这么做太笨了。更好的办法是用MathML在网页上显示公式。 MathML是XML的一个应用，是W3C推荐使用的网页数学公式显示方案。遗憾的是，目前不是所有的浏览器都支持MathML。Firefox本身包含对MathML的支持，但是你可能需要安装显示数学公式需要的字体。如果你使用Linux，你的发行版很可能包含（甚至已经安装）了该字体，请以mathml为关键词搜索以确定是否已经安装。如果你不幸使用Windows，请到这里按照说明自行安装字体。如果你更不幸使用IE浏览器，那只好自生自灭了。（开个玩笑，请安装MathPlayer插件）。 现在你可以到这里测试一下你的浏览器是不是可以正确显示MathML公式了。能显示，很好。可是，blogger并不能自动支持MathML的显示，因为blogger的网页扩展名都是html，浏览器不会自动按照xml来显示。这个问题可以这么解决：到这里下载一个小javascript文件，把这个文件直接放到你的blogger template文件的head里面，或者放到互联网上某个地方，在template的head里面加入下面一行：

<script type="text/javascript" src="javascript文件的链接"></script>

这样blogger就可以显示你的MathML公式了。 MathML语法非常复杂，要想直接用MathML写出一个公式并不容易。TtM是个把LaTeX公式转换成MathML代码的小工具。这下就方便多了。把转换好的MathML代码复制到你的blog就可以了。 下面是一个MathML显示公式的例子：

${\displaystyle U(r)=A_0\frac{W_0}{W(z)}\exp [-\frac{{\rho }^2}{W^2(z)}]\exp [-\mathrm{jkz}-\mathrm{jk}\frac{{\rho }^2}{2R(z)}+j\zeta (z)]}$

如何使用Maxima（四）

矩阵运算

1. 输入矩阵
可以用entermatrix(m,n)输入一个m乘n矩阵。如果m=n，Maxima会自动询问是什么类型的矩阵，以减少输入的工作量。

(%i1) m: entermatrix(2,2);

Is the matrix  1. Diagonal  2. Symmetric  3. Antisymmetric  4. General
Answer 1, 2, 3 or 4 :
4;
Row 1 Column 1:
a;
Row 1 Column 2:
1;
Row 2 Column 1:
0;
Row 2 Column 2:
b;

Matrix entered.
                                  [ a  1 ]
(%o1)                             [      ]
                                  [ 0  b ]

也可以用matrix以行为单位输入

(%i2) matrix([a,b,c],[1,2,3],[d,e,f]);
                                  [ a  b  c ]
                                  [         ]
(%o2)                             [ 1  2  3 ]
                                  [         ]
                                  [ d  e  f ]

2. 矩阵的转置，逆和行列式

(%i3) transpose(m);
                                   [ a  0 ]
(%o3)                              [      ]
                                   [ 1  b ]
(%i4) invert(m);
                                 [ 1     1  ]
                                 [ -  - --- ]
                                 [ a    a b ]
(%o4)                            [          ]
                                 [      1   ]
                                 [ 0    -   ]
                                 [      b   ]
(%i5) determinant(m);
(%o5)                                 a b

3. 矩阵的代数运算

(%i6) m . %o4;
                                   [ 1  0 ]
(%o6)                              [      ]
                                   [ 0  1 ]
(%i7) %o6 + %o4;
                               [ 1         1  ]
                               [ - + 1  - --- ]
                               [ a        a b ]
(%o7)                          [              ]
                               [        1     ]
                               [   0    - + 1 ]
                               [        b     ]

4. 提取矩阵的某一行或某一列

(%i8) row(m,1);
(%o8)                              [ a  1 ]
(%i9) col(%o7,2);
                                   [    1  ]
                                   [ - --- ]
                                   [   a b ]
(%o9)                              [       ]
                                   [ 1     ]
                                   [ - + 1 ]
                                   [ b     ]

5. 为矩阵增加一行（列）或几行（列）

(%i10) addcol(m,[0,1],[c,1]);
                                 [ a  1  0  c ]
(%o10)                           [            ]
                                 [ 0  b  1  1 ]
(%i11) addrow(%o10,[1,0,0,1]);
                                [ a  1  0  c ]
                                [            ]
(%o11)                          [ 0  b  1  1 ]
                                [            ]
                                [ 1  0  0  1 ]

6. 另外，+, -, *, /, ^, exp(x)等运算都可以应用于矩阵，这些运算符分别作用于每一个元素。

(%i12) m/c;
                                   [ a  1 ]
                                   [ -  - ]
                                   [ c  c ]
(%o12)                             [      ]
                                   [    b ]
                                   [ 0  - ]
                                   [    c ]

如何使用Maxima（三）

1. 复数运算
代数运算：

(%i1) a:2-3*%i$
(%i2) b:5+4*%i$
(%i3) a+b;
(%o3)                               %i + 7
(%i4) a*b;
(%o4)                        (2 - 3 %i) (4 %i + 5)

上面%i1和%i2结尾用$而不是分号;表示不输出结果。 实部和虚部：

(%i5) expand(%);
(%o5)                              22 - 7 %i
(%i6) realpart(%);
(%o6)                                 22
(%i7) imagpart(%o5);
(%o7)                                 - 7

绝对值和幅角：

(%i8) abs(%o5);
(%o8)                              sqrt(533)
(%i9) carg(%o5);
                                        7
(%o9)                             - atan(--)
                                        22

复数形式和直角形式：

(%i10) polarform(%o5);
                                     - %i atan(7/22)
(%o10)                    sqrt(533) %e
(%i11) rectform(%);
                               7                            7
(%o11)       sqrt(533) cos(atan(--)) - sqrt(533) %i sin(atan(--))
                               22                           22
(%i12) %,numer;
(%o12)                           22.0 - 7.0 %i

2. 方程数值解
找出方程f(x)=0在[a, b]范围内的数值解：

find_root(f(x), x, a, b);

如何使用Maxima（二）

1. 三角函数变换

   常用的三角函数变换有下面几种：

   trigexpand

   利用和差化积公式展开

   trigreduce

   利用积化和差等公式变成sin或cos的和

   trigsimp

   利用sin²(x) + cos² = 1等公式简化

   trigrat

   简化分数形式

   下面是一些例子：

   (%i1) sin(2*x)/cos(x)+cos(2*x);
                                 sin(2 x)
   (%o1)                         -------- + cos(2 x)
                                  cos(x)
   (%i2) trigexpand(%o1);
                                2                    2
   (%o2)                   - sin (x) + 2 sin(x) + cos (x)
   (%i3) trigreduce(%o1);
   (%o3)                     sec(x) sin(2 x) + cos(2 x)
   (%i4) trigsimp(%o1);
                             sin(2 x) + cos(x) cos(2 x)
   (%o4)                     --------------------------
                                       cos(x)
   (%i5) trigrat(%);
   (%o5)                         cos(2 x) + 2 sin(x)
   

2. 微积分

   微分：

   (%i1) f: sqrt(x-1)*exp(2*x^2);
                                                  2
                                               2 x
   (%o1)                         sqrt(x - 1) %e
   (%i2) diff(f, x);
                                                        2
                                            2        2 x
                                         2 x       %e
   (%o2)               4 sqrt(x - 1) x %e     + -------------
                                                2 sqrt(x - 1)
   (%i3) diff(sin(x)*exp(x^2), x, 2);
                          2             2                 2
                     2   x             x                 x
   (%o3)          4 x  %e   sin(x) + %e   sin(x) + 4 x %e   cos(x)
   

   不定积分：

   (%i4) integrate(%o2, x);
                                          2   log(x - 1)
                                       2 x  - ----------
                                                  2
   (%o4)                     (x - 1) %e
   

   定积分： 从这个例子可以看出，maxima是一个很聪明的代数系统，它不会像其他数学程序那样想当然的默认一个结果，而是主动问你必要的条件，然后给出正确的结果。

   (%i5) integrate(exp(-a*x^2), x, 0, inf);
   Is  a  positive, negative, or zero?
   
   positive;
                                      sqrt(%pi)
   (%o5)                              ---------
                                      2 sqrt(a)
   

   泰勒展开： 这个例子给出函数g的前三阶泰勒展开。

   (%i6) g: sin(2*x)*exp(a*x);
                                     a x
   (%o6)                           %e    sin(2 x)
   (%i7) taylor (g, x, 0, 4);
                                   2       3     3         4
                           2   (3 a  - 4) x    (a  - 4 a) x
   (%o7)/T/     2 x + 2 a x  + ------------- + ------------- + . . .
                                     3               3
   

如何使用Maxima（一）

Maxima是一个优秀的开放源代码代数运算系统，符号运算能力可以和著名商业数学软件Mathematica相媲美。这里是它的使用方法的一个简介。

1. 运行和退出

   运行：

   maxima

   进入maxima环境后，会看到开始处有(%i1)，这是行标。%i表示输入行，%o表示输出行，后面的数字表示第几个输入或输出。 退出：

   quit();

   可以看出，maxima语法和C语言类似，语句结束需要加分号。maxima里面是要区分大小写的。

2. 算术运算

   符号+, -, *, /, ^分别表示加、减、乘、除、乘方。maxima是一个精确计算的代数系统，例如分数，开方等运算会保持符号形式，不会被自动转化为小数。

   (%i1) 2/5+3/4;
                                         23
   (%o1)                                 --
                                         20
   

   如果要输出小数，可以用numer命令：

   (%i2) %,numer;
   (%o2)                                1.15
   

   百分号%代表最后一次输出的结果。另外，可以用bfloat输出任意精确度的小数，有效数字位数由变量fpprec控制，默认是16位。

   (%i3) fpprec:100;
   (%o3)                                 100
   (%i4) bfloat(%pi);
   (%o4) 3.141592653589793238462643383279502884197169399375105820974944592307816406286208998628034825342117068b0
   

   (%i3)中的冒号表示赋值。

3. 代数运算

   多项式展开：

   (%i5) (x^2+2*x*y+3*y^2)^3;
                                    2            2 3
   (%o5)                        (3 y  + 2 x y + x )
   (%i6) expand(%);
               6         5       2  4       3  3       4  2      5      6
   (%o6)   27 y  + 54 x y  + 63 x  y  + 44 x  y  + 21 x  y  + 6 x  y + x
   

   等量代换：

   (%i7) %o6, x=2/z;
                    5        4        3        2
               108 y    252 y    352 y    336 y    192 y   64       6
   (%o7)       ------ + ------ + ------ + ------ + ----- + -- + 27 y
                 z         2        3        4       5      6
                          z        z        z       z      z
   

   提取公分母：

   (%i8) ratsimp(%);
             6  6        5  5        4  4        3  3        2  2
         27 y  z  + 108 y  z  + 252 y  z  + 352 y  z  + 336 y  z  + 192 y z + 64
   (%o8) -----------------------------------------------------------------------
                                            6
                                           z
   

   因式分解：

   (%i9) factor(%);
                                   2  2             3
                               (3 y  z  + 4 y z + 4)
   (%o9)                       ----------------------
                                          6
                                         z
   

   解非线性方程组：

   (%i10) a+b+2*c=0;
   (%o10)                          2 c + b + a = 0
   (%i11) 2*a-b=1;
   (%o11)                            2 a - b = 1
   (%i12) a-2*b*c=5;
   (%o12)                           a - 2 b c = 5
   (%i13) solve([%o10,%o11,%o12],[a,b,c]);
                sqrt(7) + 1      2 sqrt(7) - 1      sqrt(7) - 14
   (%o13) [[a = -----------, b = -------------, c = -------------],
                     3                 3            4 sqrt(7) - 2
                             sqrt(7) - 1        2 sqrt(7) + 1      sqrt(7) + 14
                      [a = - -----------, b = - -------------, c = -------------]]
                                  3                   3            4 sqrt(7) + 2
   

Monte Carlo方法

The Review of Particle Physics有一篇关于Monte Carlo方法的Review，在Mathematical Tools里面，介绍了Monte Carlo方法和任意分布随机数的生成，地址是 http://pdg.lbl.gov/2005/reviews/monterpp.pdf。 另外还有一些关于概率和统计的Review，对进行计算机模拟实验很有帮助。