Matlab数值计算

数值计算

函数句柄

在 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^{T}A$可逆，则给出最小二乘解。

>> A = [2, 1; 1, 3];
>> b = [8; 13];
>> x = A \ bx =2.20003.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.8955fval =-2.2204e-16exitflag =1output = 包含以下字段的 struct:intervaliterations: 11iterations: 6funcCount: 29algorithm: '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.7071fval =1.0e-11 *0.22820exitflag =1output = 包含以下字段的 struct:iterations: 4funcCount: 15algorithm: 'trust-region-dogleg'firstorderopt: 3.2275e-12message: '方程已解。...'jacobian =1.4142    1.41421.0000   -1.0000

4. roots

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

roots(c)

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

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

函数极值的求解

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 =2fval =0exitflag =1output = 包含以下字段的 struct:iterations: 5funcCount: 6algorithm: '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⩽b$
• Aeq, beq线性等式约束 $Aeqx=beq$
• lb, ub变量下界和上界
• nonlcon非线性约束函数，返回 [c(x), ceq(x)]
• options优化参数

如果某些约束没有，可用 [] 占位
示例：
最小化目标函数：
$$\begin{array}{c}f(x_1,x_2)=(x_1-2{)}^2+(x_2-1{)}^2\\ s.t.\{\begin{array}{l}{x}_1+x_2⩽3\\ x_1-x_2=0\\ x_1^2+x_2^2-4⩽0\\ \sin (x_1+x_2)-1=0\\ 0⩽x_1⩽2,0⩽x_2⩽2\end{array}\end{array}$$

% 非线性约束函数文件nonlcon.m
function [c, ceq] = nonlcon(x)% 非线性不等式：c(x) <= 0c = x(1)^2 + x(2)^2 - 4;% 非线性等式：ceq(x) = 0ceq = 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.78540.7854fval =1.5213exitflag =2output = 包含以下字段的 struct:iterations: 27funcCount: 93constrviolation: 2.2204e-16stepsize: 8.0922e-11algorithm: 'interior-point'firstorderopt: 0.0814cgiterations: 0message: '可能存在局部最小值。满足约束。...'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.0000fval =-0.1016exitflag =1output = 包含以下字段的 struct:iterations: 78funcCount: 148algorithm: 'Nelder-Mead simplex direct search'message: '优化已终止:...'>> [x, fval, exitflag, output] = fminunc(fun, x0)找到局部最小值。优化已完成，因为梯度大小小于
最优性容差的值。<停止条件详细信息>x =3.4728   -1.0000fval =-0.1016exitflag =1output = 包含以下字段的 struct:iterations: 7funcCount: 24stepsize: 6.9597e-05lssteplength: 1firstorderopt: 2.3117e-07algorithm: '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.5584err =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     35     5>> dA = diff(A, 1, 2)dA =111

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 =-3den =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.6400S = 包含以下字段的 struct:R: [2x2 double]df: 3normr: 0.1643rsquared: 0.9952mu =3.00001.5811