草庐IT

【最优化理论】牛顿法+Matlab代码实现

WSKH0929 2023-08-30 原文

文章目录


1 牛顿法简介

牛顿迭代法(Newton’s method)又称为牛顿-拉夫逊(拉弗森)方法(Newton-Raphson method),它是牛顿在17世纪提出的一种在实数域和复数域上近似求解方程的方法。

多数方程不存在求根公式,因此求精确根非常困难,甚至不可解,从而寻找方程的近似根就显得特别重要。方法使用函数 f ( x ) f(x) f(x) 的泰勒级数的前面几项来寻找方程 f ( x ) = 0 f(x)=0 f(x)=0 的根。牛顿迭代法是求方程根的重要方法之一,其最大优点是在方程 f ( x ) = 0 f(x)=0 f(x)=0 的单根附近具有平方收敛,而且该法还可以用来求方程的重根、复根,此时线性收敛,但是可通过一些方法变成超线 性收㪉。另外该方法广泛用于计算机编程中。


2 牛顿法原理

r r r f ( x ) = 0 f(x)=0 f(x)=0 的根,选取 x 0 x_{0} x0 作为 r r r 的初始近似值,过点 ( x 0 , f ( x 0 ) ) \left(x_{0}, f\left(x_{0}\right)\right) (x0,f(x0)) 做曲线 y = f ( x ) y=f(x) y=f(x) 的切线 L L L
L : y = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) L: y=f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right) L:y=f(x0)+f(x0)(xx0) ,则 L L L x x x 轴交点的横坐标 x 1 = x 0 − f ( x 0 ) f ′ ( x 0 ) x_{1}=x_{0}-\frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)} x1=x0f(x0)f(x0) ,称 x 1 x_{1} x1 r r r 的一次近似值。过点 ( x 1 , f ( x 1 ) ) \left(x_{1}, f\left(x_{1}\right)\right) (x1,f(x1)) 做曲线 y = f ( x ) y=f(x) y=f(x) 的切线,并求该切线与 × \times × 轴交点的横坐标 x 2 = x 1 − f ( x 1 ) f ′ ( x 1 ) x_{2}=x_{1}-\frac{f\left(x_{1}\right)}{f^{\prime}\left(x_{1}\right)} x2=x1f(x1)f(x1) ,称 x 2 x_{2} x2 r \mathrm{r} r 的二次近似值。重曷 以上过程,得 r r r 的近似值序列,其中, x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1}=x_{n}-\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)} xn+1=xnf(xn)f(xn) 称为 r r r n + 1 n+1 n+1 次近似值,上式称为牛顿迭代公式。

用牛顿迭代法解非线性方程,是把非线性方程 f ( x ) = 0 f(x)=0 f(x)=0 线性化的一种近似方法。把 f ( x ) f(x) f(x) 在点 x 0 x_{0} x0 的桌邻域内展开成泰勒 级数 f ( x ) = f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) + f ′ ′ ( x 0 ) ( x − x 0 ) 2 2 ! + ⋯ + f ( n ) ( x 0 ) ( x − x 0 ) n n ! + R n ( x ) f(x)=f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right)+\frac{f^{\prime \prime}\left(x_{0}\right)\left(x-x_{0}\right)^{2}}{2 !}+\cdots+\frac{f^{(n)}\left(x_{0}\right)\left(x-x_{0}\right)^{n}}{n !}+R_{n}(x) f(x)=f(x0)+f(x0)(xx0)+2!f′′(x0)(xx0)2++n!f(n)(x0)(xx0)n+Rn(x) ,取其线性部分 (即泰勒展开的前两项),并令其等于 0 ,即 f ( x 0 ) + f ′ ( x 0 ) ( x − x 0 ) = 0 f\left(x_{0}\right)+f^{\prime}\left(x_{0}\right)\left(x-x_{0}\right)=0 f(x0)+f(x0)(xx0)=0 ,以此作为非线性方程 f ( x ) = 0 f(x)=0 f(x)=0 的近似方程, 若 f ′ ( x 0 ) ≠ 0 f^{\prime}\left(x_{0}\right) \neq 0 f(x0)=0 ,则其解为 x 1 = x 0 − f ( x 0 ) f ′ ( x 0 ) x_{1}=x_{0}-\frac{f\left(x_{0}\right)}{f^{\prime}\left(x_{0}\right)} x1=x0f(x0)f(x0) ,这样,得到牛顿迭代法的一个朱代关系式: x n + 1 = x n − f ( x n ) f ′ ( x n ) x_{n+1}=x_{n}-\frac{f\left(x_{n}\right)}{f^{\prime}\left(x_{n}\right)} xn+1=xnf(xn)f(xn)

已经证明,如果是连续的,并且待求的零点是孤立的,那么在零点周围存在一个区域,只要初始值位于这个邻近区域内,那 么牛顿法必定收敛。并且,如果不为 0 , 那么牛顿法将具有平方收敛的性能. 粗略的说,这意味着每造代一次,牛顿法结果的有效 数字将增加一倍。


3 牛顿法推导




4 Matlab代码实现

下面用Matlab代码求解上面的示例。

clear;clc;

% 定义原函数
syms xx yy
fy(xx,yy) = 0.5 * xx^2 + 2 * yy^2;

% 确定迭代次数
n = 10
% 确定初始点
x0 = 1
y0 = 1
% 求初始点函数值
fy(x0,y0)
% 求函数梯度
xf = -5:0.2:5;
yf = xf';
ff = 0.5 * xf.^2 + 2 * yf.^2;
surf(xf,yf,ff)
xlabel('x')
ylabel('y')
zlabel('z')
view([119.1 40.8])
[fx,fy] = gradient(ff,0.2);

% 提取点初始点处的梯度值
t = (xf == x0) & (yf == y0);
indt = find(t);
f_grad = [fx(indt) fy(indt)]
% 求海森矩阵
syms x y
f(x,y) = 0.5 * x^2 + 2 * y^2;
H = hessian(f,[x,y])
% 迭代
for i=1:n

    % 判断是否可以跳出(如果梯度向量都接近0,就跳出)
    b = 0;
    for j = 1:length(f_grad)
        if f_grad(j) > 0.000001
            b = 1;
            break
        end
    end
    if b==0
        break
    end

    % 确定下降方向
    d = -inv(H)*(f_grad)';
    dk = d(x0,y0);
    
    % 确定步长,牛顿法步长为1
    a = 1;

    % 获取下一状态的点
    newX = [x0,y0] + dk' .* a
    x0 = newX(1);
    y0 = newX(2);

    % 更新梯度信息
    t = (xf == x0) & (yf == y0);
    indt = find(t);
    f_grad = [fx(indt) fy(indt)];

end


5 低版本Matlab报错

最近有朋友向我反应代码运行会报错,具体报错内容如下:

他使用的matlab版本是2016a,推测可能是低版本不支持(151)的矩阵和(511)的矩阵直接做运算,如果大家有遇到这样的报错的话,可以试一下将原代码的16、17行删去,换成以下代码应该就可以了:

yf = xf;
s = size(xf,2);
ff = zeros(s,s);
for i = 1 : s
    for j = 1 : s
        obj = 0.5 * xf(i)^2 + 2*yf(j)^2;
        ff(i,j) = obj;
    end
end

有关【最优化理论】牛顿法+Matlab代码实现的更多相关文章

  1. ruby - 如何在 buildr 项目中使用 Ruby 代码? - 2

    如何在buildr项目中使用Ruby?我在很多不同的项目中使用过Ruby、JRuby、Java和Clojure。我目前正在使用我的标准Ruby开发一个模拟应用程序,我想尝试使用Clojure后端(我确实喜欢功能代码)以及JRubygui和测试套件。我还可以看到在未来的不同项目中使用Scala作为后端。我想我要为我的项目尝试一下buildr(http://buildr.apache.org/),但我注意到buildr似乎没有设置为在项目中使用JRuby代码本身!这看起来有点傻,因为该工具旨在统一通用的JVM语言并且是在ruby中构建的。除了将输出的jar包含在一个独特的、仅限ruby​​

  2. ruby-on-rails - Rails 源代码 : initialize hash in a weird way? - 2

    在rails源中:https://github.com/rails/rails/blob/master/activesupport/lib/active_support/lazy_load_hooks.rb可以看到以下内容@load_hooks=Hash.new{|h,k|h[k]=[]}在IRB中,它只是初始化一个空哈希。和做有什么区别@load_hooks=Hash.new 最佳答案 查看rubydocumentationforHashnew→new_hashclicktotogglesourcenew(obj)→new_has

  3. ruby - 如何根据特征实现 FactoryGirl 的条件行为 - 2

    我有一个用户工厂。我希望默认情况下确认用户。但是鉴于unconfirmed特征,我不希望它们被确认。虽然我有一个基于实现细节而不是抽象的工作实现,但我想知道如何正确地做到这一点。factory:userdoafter(:create)do|user,evaluator|#unwantedimplementationdetailshereunlessFactoryGirl.factories[:user].defined_traits.map(&:name).include?(:unconfirmed)user.confirm!endendtrait:unconfirmeddoenden

  4. ruby-on-rails - 浏览 Ruby 源代码 - 2

    我的主要目标是能够完全理解我正在使用的库/gem。我尝试在Github上从头到尾阅读源代码,但这真的很难。我认为更有趣、更温和的踏脚石就是在使用时阅读每个库/gem方法的源代码。例如,我想知道RubyonRails中的redirect_to方法是如何工作的:如何查找redirect_to方法的源代码?我知道在pry中我可以执行类似show-methodmethod的操作,但我如何才能对Rails框架中的方法执行此操作?您对我如何更好地理解Gem及其API有什么建议吗?仅仅阅读源代码似乎真的很难,尤其是对于框架。谢谢! 最佳答案 Ru

  5. ruby - 模块嵌套代码风格偏好 - 2

    我的假设是moduleAmoduleBendend和moduleA::Bend是一样的。我能够从thisblog找到解决方案,thisSOthread和andthisSOthread.为什么以及什么时候应该更喜欢紧凑语法A::B而不是另一个,因为它显然有一个缺点?我有一种直觉,它可能与性能有关,因为在更多命名空间中查找常量需要更多计算。但是我无法通过对普通类进行基准测试来验证这一点。 最佳答案 这两种写作方法经常被混淆。首先要说的是,据我所知,没有可衡量的性能差异。(在下面的书面示例中不断查找)最明显的区别,可能也是最著名的,是你的

  6. ruby - 寻找通过阅读代码确定编程语言的ruby gem? - 2

    几个月前,我读了一篇关于ruby​​gem的博客文章,它可以通过阅读代码本身来确定编程语言。对于我的生活,我不记得博客或gem的名称。谷歌搜索“ruby编程语言猜测”及其变体也无济于事。有人碰巧知道相关gem的名称吗? 最佳答案 是这个吗:http://github.com/chrislo/sourceclassifier/tree/master 关于ruby-寻找通过阅读代码确定编程语言的rubygem?,我们在StackOverflow上找到一个类似的问题:

  7. ruby - Net::HTTP 获取源代码和状态 - 2

    我目前正在使用以下方法获取页面的源代码:Net::HTTP.get(URI.parse(page.url))我还想获取HTTP状态,而无需发出第二个请求。有没有办法用另一种方法做到这一点?我一直在查看文档,但似乎找不到我要找的东西。 最佳答案 在我看来,除非您需要一些真正的低级访问或控制,否则最好使用Ruby的内置Open::URI模块:require'open-uri'io=open('http://www.example.org/')#=>#body=io.read[0,50]#=>"["200","OK"]io.base_ur

  8. 程序员如何提高代码能力? - 2

    前言作为一名程序员,自己的本质工作就是做程序开发,那么程序开发的时候最直接的体现就是代码,检验一个程序员技术水平的一个核心环节就是开发时候的代码能力。众所周知,程序开发的水平提升是一个循序渐进的过程,每一位程序员都是从“菜鸟”变成“大神”的,所以程序员在程序开发过程中的代码能力也是根据平时开发中的业务实践来积累和提升的。提高代码能力核心要素程序员要想提高自身代码能力,尤其是新晋程序员的代码能力有很大的提升空间的时候,需要针对性的去提高自己的代码能力。提高代码能力其实有几个比较关键的点,只要把握住这些方面,就能很好的、快速的提高自己的一部分代码能力。1、多去阅读开源项目,如有机会可以亲自参与开源

  9. Matlab imread()读到了什么 (浅显 当复习文档了) - 2

    matlab打开matlab,用最简单的imread方法读取一个图像clcclearimg_h=imread('hua.jpg');返回一个数组(矩阵),往往是a*b*cunit8类型解释一下这个三维数组的意思,行数、数和层数,unit8:指数据类型,无符号八位整形,可理解为0~2^8的数三个层数分别代表RGB三个通道图像rgb最常用的是24-位实现方法,即RGB每个通道有256色阶(2^8)。基于这样的24-位RGB模型的色彩空间可以表现256×256×256≈1670万色当imshow传入了一个二维数组,它将以灰度方式绘制;可以把图像拆分为rgb三层,可以以灰度的方式观察它figure(1

  10. 华为OD机试用Python实现 -【明明的随机数】 2023Q1A - 2

    华为OD机试题本篇题目:明明的随机数题目输入描述输出描述:示例1输入输出说明代码编写思路最近更新的博客华为od2023|什么是华为od,od薪资待遇,od机试题清单华为OD机试真题大全,用Python解华为机试题|机试宝典【华为OD机试】全流程解析+经验分享,题型分享,防作弊指南华为o

随机推荐