数值计算
函数句柄
在 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 A A可逆,则x=A \ b
等价于x=inv(A)*b
,若 A A A的行数大于列数且 A T A A^TA ATA可逆,则给出最小二乘解。
>> 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 f(a)f(b) < 0 f(a)f(b)<0,函数在区间中必须变号)
options
是求解选项
输出:
x
是求得的近似解fval
是函数fun
在x
处的取值exitflag
表示退出原因:- 1 1 1表示找到解
- − 1 -1 −1表示未收敛
- − 3 -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 ] [x_1, x_2] [x1,x2]上的最小值点与最小值。
示例:
>> 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
线性不等式约束 A x ⩽ b Ax\leqslant b Ax⩽bAeq, beq
线性等式约束 A e q x = b e q Aeq x = beq Aeqx=beqlb, ub
变量下界和上界nonlcon
非线性约束函数,返回[c(x), ceq(x)]
options
优化参数
如果某些约束没有,可用 []
占位
示例:
最小化目标函数:
f ( x 1 , x 2 ) = ( x 1 − 2 ) 2 + ( x 2 − 1 ) 2 s.t. { 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 \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*} f(x1,x2)=(x1−2)2+(x2−1)2s.t.⎩
⎨
⎧x1+x2⩽3x1−x2=0x12+x22−4⩽0sin(x1+x2)−1=00⩽x1⩽2,0⩽x2⩽2
% 非线性约束函数文件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 10^{-6} 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