Matlab数值计算

MATLAB数值计算

数值计算
- 函数句柄
- 匿名函数
- 线性与非线性方程组求解
- - 1. `\`（左除运算）
  - 2. `fzero`
  - 3. `fsolve`
  - 4. `roots`
- 函数极值的求解
- - 1. `fminbnd`
  - 2. `fmincon`
  - 3. `fminsearch`与`fminunc`
- 数值积分
- - 1. `quad` / `quadl`
  - 2. `quadgk`
  - 3. `integral`
  - 4. `trapz`
  - 5. `dblquad`, `quad2d`, `integral2`
  - 6. `triplequad`, `integral3`
- 数值微分
- - 1. `diff`
  - 2. `polyder`
  - 3. `polyval`
  - 4. `spline`
  - 5. `polyfit`

数值计算

函数句柄

在 MATLAB 中，当你用文件定义一个函数，比如：

function f = myfun(x)
    f = sin(x) - x/2;
end

这表示你创建了一个用户自定义函数 myfun，输入为 x，输出为 f。你需要将这段代码保存为文件，文件名必须是 myfun.m，并确保该文件在 MATLAB 的当前路径或工作路径下。
如何调用这个函数？
方法一：直接传值调用

myfun(1.5)

ans =

    0.2475

方法二：使用@

x = fzero(@myfun, 1);

@myfun 是一个函数句柄，告诉 MATLAB 使用myfun函数。

匿名函数

语法：

@(x) f(x)

@后面括号中是自变量，f(x)处写函数具体形式，结果表示一个函数。若f(x)处写的是函数向量，则结果表示函数组，如fun = @(x) [x(1)^2 + x(2)^2 - 1; x(1) - x(2)]

线性与非线性方程组求解

1. `\`（左除运算）

用于解线性方程组 Ax = b。若矩阵 $A$ 可逆，则x=A \ b等价于x=inv(A)*b，若 $A$ 的行数大于列数且 $A^TA$ 可逆，则给出最小二乘解。

>> A = [2, 1; 1, 3];
>> b = [8; 13];
>> x = A \ b

x =

    2.2000
    3.6000

2. `fzero`

求解非线性单变量方程 f(x) = 0。
语法：

[x, fval, exitflag, output] = fzero(fun, x0, options)

输入：
fun 是目标函数，可以是：

字符串句柄，例如：'sin(x)-x/2'
函数句柄，例如：@(x) sin(x) - x/2

x0 是初始猜测值，可以是：

单个数（算法从该点出发）
一个区间 [a, b]（要求 $f (a) f (b) < 0$ ，函数在区间中必须变号）

options是求解选项

输出：

x是求得的近似解
fval是函数fun在x处的取值
exitflag表示退出原因：
1. $1$ 表示找到解
2. $- 1$ 表示未收敛
3. $- 3$ 表示函数在区间上不变号
output是结构体，包含迭代信息，如迭代次数、函数调用次数等

示例：

>> fun = @(x) sin(x) - x/2;
>> x0 = 1;
>> [x, fval, exitflag, output] = fzero(fun, x0)

x =

    1.8955


fval =

  -2.2204e-16


exitflag =

     1


output = 

  包含以下字段的 struct:

    intervaliterations: 11
            iterations: 6
             funcCount: 29
             algorithm: 'bisection, interpolation'
               message: '在区间 [0.0949033, 1.9051] 中发现零'

3. `fsolve`

求解非线性方程组。
语法：

[x, fval, exitflag, output, jacobian] = fsolve(fun, x0, options)

jacobian表示在解处的Jacobian矩阵

>> fun = @(x) [x(1)^2 + x(2)^2 - 1; x(1) - x(2)];
>> x0 = [0.5, 0.5];
>> [x, fval, exitflag, output, jacobian] = fsolve(fun, x0)

方程已解。

fsolve 已完成，因为按照函数容差的值衡量，
函数值向量接近于零，并且按照梯度的值衡量，
问题似乎为正则问题。

<停止条件详细信息>

x =

    0.7071    0.7071


fval =

   1.0e-11 *

    0.2282
         0


exitflag =

     1


output = 

  包含以下字段的 struct:

       iterations: 4
        funcCount: 15
        algorithm: 'trust-region-dogleg'
    firstorderopt: 3.2275e-12
          message: '方程已解。...'


jacobian =

    1.4142    1.4142
    1.0000   -1.0000

4. `roots`

求解一元多项式的所有根。
语法：

roots(c)

c是多项式的系数向量，按降幂排序

>> c = [1 -3 2];  % 表示 x^2 - 3x + 2
>> r = roots(c)

r =

     2
     1

>>

函数极值的求解

MATLAB只求最小值

1. `fminbnd`

求解有界有约束单变量函数的最小值，约束体现在是求解区间上的最小值。
语法：

[x, fval, exitflag, output]=fminbnd(fun, x1, x2, options)

给出fun在 $x_1, x_2]$ 上的最小值点与最小值。
示例：

>> f = @(x) (x-2).^2;
>> [x, fval, exitflag, output] = fminbnd(f, 0, 4)

x =

     2


fval =

     0


exitflag =

     1


output = 

  包含以下字段的 struct:

    iterations: 5
     funcCount: 6
     algorithm: 'golden section search, parabolic interpolation'
       message: '优化已终止:...'

2. `fmincon`

求解有约束多变量函数的最小值。
语法：

[x,fval] = fmincon(fun, x0, A, b, Aeq, beq, lb, ub, nonlcon, options)

fun目标函数
x0初始点（必须给）
A, b线性不等式约束 $Ax\leqslant b$
Aeq, beq线性等式约束 $A e q x = b e q$
lb, ub变量下界和上界
nonlcon非线性约束函数，返回 [c(x), ceq(x)]
options优化参数

如果某些约束没有，可用 [] 占位
示例：
最小化目标函数：
$\begin{gather*} f(x_1,x_2)=(x_1-2)^2+(x_2-1)^2 \\ \operatorname{s.t.} \begin{cases} x_1+x_2\leqslant3 \\ x_1-x_2=0 \\ x_1^2+x_2^2-4\leqslant0 \\ \sin(x_1+x_2)-1=0 \\ 0\leqslant x_1\leqslant2,\quad0\leqslant x_2\leqslant2 \end{cases} \end{gather*}$

% 非线性约束函数文件nonlcon.m
function [c, ceq] = nonlcon(x)
    % 非线性不等式：c(x) <= 0
    c = x(1)^2 + x(2)^2 - 4;
    % 非线性等式：ceq(x) = 0
    ceq = sin(x(1) + x(2)) - 1;
end

% 目标函数
>> fun = @(x) (x(1) - 2)^2 + (x(2) - 1)^2;

% 初始点
>> x0 = [1; 1];

% 线性不等式约束 A*x <= b
>> A = [1, 1];
>> b = 3;

% 线性等式约束 Aeq*x = beq
>> Aeq = [1, -1];
>> beq = 0;

% 变量上下界
>> lb = [0; 0];
>> ub = [2; 2];

% 求解
>> [x_opt, fval, exitflag, output] = fmincon(fun, x0, A, b, Aeq, beq, lb, ub, @nonlcon)

找到具有较低目标函数值的可行点，但不满足最优性条件。请参阅 output.bestfeasible。


可能存在局部最小值。满足约束。

fmincon 已停止，因为当前步长小于
步长容差值并且在约束容差值范围内满足约束。

<停止条件详细信息>

x_opt =

    0.7854
    0.7854


fval =

    1.5213


exitflag =

     2


output = 

  包含以下字段的 struct:

         iterations: 27
          funcCount: 93
    constrviolation: 2.2204e-16
           stepsize: 8.0922e-11
          algorithm: 'interior-point'
      firstorderopt: 0.0814
       cgiterations: 0
            message: '可能存在局部最小值。满足约束。...'
       bestfeasible: [1x1 struct]

3. `fminsearch`与`fminunc`

二者是求无约束多变量最小值问题的函数，其中fminsearch使用Nelder-Mead法。

[x,fval,exitflag,output] = fminsearch(fun, x0, options)
[x,fval,exitflag,output] = fminunc(fun, x0, options)

>> fun = @(x) (x(1) - 3)^2 + (x(2) + 1)^2 + sin(x(1));
>> x0 = [0, 0];
>> [x, fval, exitflag, output] = fminsearch(fun, x0)

x =

    3.4728   -1.0000


fval =

   -0.1016


exitflag =

     1


output = 

  包含以下字段的 struct:

    iterations: 78
     funcCount: 148
     algorithm: 'Nelder-Mead simplex direct search'
       message: '优化已终止:...'

>> [x, fval, exitflag, output] = fminunc(fun, x0)

找到局部最小值。

优化已完成，因为梯度大小小于
最优性容差的值。

<停止条件详细信息>

x =

    3.4728   -1.0000


fval =

   -0.1016


exitflag =

     1


output = 

  包含以下字段的 struct:

       iterations: 7
        funcCount: 24
         stepsize: 6.9597e-05
     lssteplength: 1
    firstorderopt: 2.3117e-07
        algorithm: 'quasi-newton'
          message: '找到局部最小值。...'

数值积分

1. `quad` / `quadl`

自适应辛普森法与自适应 Lobatto 高阶法求积分
语法：

[q, n] = quad(fun, a, b, tol)
[q, n] = quadl(fun, a, b)

q返回积分值，n返回函数计算的次数，tol表示精度，默认为 $10^{-6}$ ，上式表示计算函数fun从a到b的积分
示例：

>> f = @(x) x.^2;
>> quad(f, 0, 1)

ans =

    0.3333

>> quadl(f, 0, 1)

ans =

    0.3333

2. `quadgk`

高精度 Gauss-Lobatto 法
语法

[q, err] = quadgk(fun, a, b)

err表示误差估计
示例：

>> f = @(x) sin(x)./x;
>> [q, err] = quadgk(f, 0.1, 10)

q =

    1.5584


err =

   1.7052e-16

3. `integral`

全局自适应辛普森公式
语法：

q = integral(fun, a, b)

示例：

>> f = @(x) exp(-x.^2);
>> I = integral(f, -Inf, Inf)

I =

    1.7725

4. `trapz`

梯形公式
语法：

I = trapz(x, y)

示例：

>> x = 0:0.1:10;
>> y = sin(x);
>> I = trapz(x, y)

I =

    1.8375

5. `dblquad`, `quad2d`, `integral2`

计算二重积分。
语法：

q = dblquad(fun, xmin, xmax, ymin, ymax, tol, method)
q = quad2d(fun, xmin, xmax, ymin, ymax)
q = integral2(fun, xmin, xmax, ymin, ymax)

示例：

>> f = @(x, y) x .* y;
>> q = dblquad(f, 0, 1, 0, 1)

q =

    0.2500

>> q = quad2d(f, 0, 1, 0, 1) 

q =

    0.2500

>> q = integral2(f, 0, 1, 0, 1)

q =

    0.2500

6. `triplequad`, `integral3`

三重积分。
语法：

q = triplequad(fun, xmin, xmax, ymin, ymax, zmin, zmax)
q = integral3(fun, xmin, xmax, ymin, ymax, zmin, zmax)

示例：

>> f = @(x, y, z) x + y + z;
>> q = integral3(f, 0, 1, 0, 1, 0, 1)

q =

    1.5000

>> q = triplequad(f, 0, 1, 0, 1, 0, 1)

q =

    1.5000

数值微分

1. `diff`

向前差分。
语法：

% 计算向量 x 的一阶向前差分：x(i+1) - x(i)
diff(x, n)
% 对矩阵 A 的第 dim 个维度计算 n 阶向前差分，dim = 1（默认）：按列方向差分，dim = 2：按行方向差分
diff(A, n, dim)

示例：

>> x = [1 2 4 7];
>> dx = diff(x)

dx =

     1     2     3

>> dx = diff(x, 2)

dx =

     1     1

>> A = [1 2; 4 5; 9 10];
>> dA = diff(A, 1, 1)

dA =

     3     3
     5     5

>> dA = diff(A, 1, 2)

dA =

     1
     1
     1

2. `polyder`

求多项式导函数。

语法：

dp = polyder(p)                 % 求多项式 p 的导函数
dp = polyder(p, q)              % 求 p*q 的导函数
[dp, dq] = polyder(p, q)        % 求 p/q 的导函数，dp返回分子多项式的系数，dq返回分母多项式的系数

示例：

>> p = [3 0 -4];         % 表示 3x^2 - 4
>> dp = polyder(p)       % 导数为 6x 

dp =

     6     0

>> p = [1 2];            % p(x) = x + 2
>> q = [3 0 -1];         % q(x) = 3x^2 - 1
>> dpq = polyder(p, q)   % 求 (p*q)的导函数

dpq =

     9    12    -1

>> p = [1 2];                % 分子 p(x) = x + 2
>> q = [1 -1];               % 分母 q(x) = x - 1
>> [num, den] = polyder(p, q)

num =

    -3


den =

     1    -2     1

3. `polyval`

多项式求值。

p = [1 -3 2];
y = polyval(p, 2)

y =

     0

4. `spline`

三次样条插值。
语法：

yi = spline(x, y, xi)

x, y 是已知数据点
xi 是需要插值的位置
返回 yi = f(xi)

示例：

>> x = 0:10;
>> y = sin(x);
>> xi = 0:0.1:1;
>> yi = spline(x, y, xi)

yi =

         0    0.1118    0.2181    0.3187    0.4134    0.5017    0.5837    0.6588    0.7270    0.7880    0.8415

5. `polyfit`

数据拟合为多项式。
语法：

[p, S, mu] = polyfit(x, y, n)

拟合n次多项式，返回误差结构体，mu中第一个元素为x的均值，第二个元素为x的标准差
示例：

>> x = 1:5;
>> y = [2.2 2.8 3.6 4.5 5.1];
>> [p, S, mu] = polyfit(x, y, 1)

p =

    1.1859    3.6400


S = 

  包含以下字段的 struct:

           R: [2x2 double]
          df: 3
       normr: 0.1643
    rsquared: 0.9952


mu =

    3.0000
    1.5811