ICode9

精准搜索请尝试: 精确搜索
首页 > 编程语言> 文章详细

matlab练习程序(数值雅克比)

2022-02-02 13:33:00  阅读:175  来源: 互联网

标签:pre 雅克 epsilon inv 练习 数值 matlab delta


数值雅克比本质就是对函数的每一维分别做数值微分,再组合为雅克比矩阵即可。

通常我们做最优化的时候要计算函数的雅克比矩阵,但是如果函数的解析式比较复杂,求其偏导解析解会非常麻烦。

虽然可以利用Mathematica或者Matlab的符号运算进行求解,不过有时候得到的解析解也是很复杂的,再转写成代码如果错一个符号很可能就找不到错误来源了。

利用数值方法求偏导,得到雅克比矩阵,在进行优化求解,通常可以简化整个过程。

下面用之前文章中高斯牛顿法优化作为例子,来实践一下如何计算数值雅克比。

matlab代码如下:

clear all;close all;clc;

epsilon = 1e-6;
epsilon_inv = 1/epsilon;

a=1;b=2;c=1;                        %待求解的系数
x = (0:0.01:1)';
w = rand(length(x),1)*2-1;          %生成噪声
y = exp(a*x.^2+b*x+c) + w;          %带噪声的模型 
plot(x,y,'.')

pre=rand(3,1);                      %步骤1
for i=1:100
    
    f0 = exp(pre(1)*x.^2+pre(2)*x+pre(3));
    g = y - f0;                     %步骤2中的误差 
         
    %数值计算雅克比
    for k = 1:length(pre)
        x_ = pre;
        x_(k) = pre(k) + epsilon;
        jac(:, k) = (exp(x_(1)*x.^2+x_(2)*x+x_(3)) - f0) .* epsilon_inv;
    end
    J = jac;
    
    delta = inv(J'*J)*J'* g;        %步骤3,inv(J'*J)*J'为H的逆
    
    pcur = pre+delta;               %步骤4
    if norm(delta) <1e-16
        break;
    end
    pre = pcur;
end

hold on;
plot(x,exp(a*x.^2+b*x+c),'r');
plot(x,exp(pre(1)*x.^2+pre(2)*x+pre(3)),'g');

%比较一下
[a b c]
pre'

结果如下:

标签:pre,雅克,epsilon,inv,练习,数值,matlab,delta
来源: https://www.cnblogs.com/tiandsp/p/15860261.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有