Matlab数值计算

发布于:2025-06-04 ⋅ 阅读:(23) ⋅ 点赞:(0)

数值计算

函数句柄

在 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是函数funx处的取值
  • exitflag表示退出原因:
    1. 1 1 1表示找到解
    2. − 1 -1 1表示未收敛
    3. − 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 Axb
  • Aeq, beq线性等式约束 A e q x = b e q Aeq x = beq Aeqx=beq
  • lb, 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)=(x12)2+(x21)2s.t. x1+x23x1x2=0x12+x2240sin(x1+x2)1=00x12,0x22

% 非线性约束函数文件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. fminsearchfminunc

二者是求无约束多变量最小值问题的函数,其中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} 106,上式表示计算函数funab的积分
示例:

>> 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
    

网站公告

今日签到

点亮在社区的每一天
去签到