Matlab 非线性迭代法（1）牛顿法

一、牛顿迭代公式

1、定义

 2、原理推导

 泰勒公式：

 常用的8个泰勒公式：

 推导：

将f(x)=f(x)在Xk 处的泰勒公式展开：

f(x)=f(Xk)+f(Xk)`(X-Xk)+    f(Xk)`/2  *(x-Xk)^2+..........

我们吧线性的一部分先拿出来：f(x)=f(Xk)+f(Xk)`(X-Xk) =0；

 3、几何意义

盗图：

MATLAB之牛顿迭代法_newton迭代法 matlab_天涯铭的博客-CSDN博客

 图中直线表示的是在x0出的斜率，可以表示为：

f(x)`=f(x)/x0-x1

化简后得到：  x1=x0-f(x0)/f(x)`

4、牛顿法的收敛条件

若 x在【a，b】区间是2连续的，那么牛顿法的条件如下：

a、f(a)*f(b)<0;

b、f(x)`!=0

c、f(x)``在【a,b】区间上是同符号

d、f(x)``*f(x0)>0  .x0在【a,b】;

牛顿迭代法的收敛性取决于x0的选择，所以x0的选取很重要

 

 二、算法流程

 eg: 计算

9.*x.^2-sin(x)-1

的极值点 

先画出这个曲线

那我们做的时候可以将初始化设定为2

导数：

clc

syms x
f=9*x^2-sin(x)-1;
jacobian(f,x)

% f 表示的是函数 
% x_init 表示的是我们 猜想的一个值  
%  epsilon 表示的 是收敛精度
%  nummaxitator 最大的迭代次数
function [x_ok,y_ok,num_itorator]=Albert_NewTon_Method(f,x_init,epsilon,nummaxitator)
format long;
if nargin==2
    epsilon=0.000001;
    nummaxitator=200;
end

%  计算函数f 的一阶导数
first_f=diff(f);
%  计算函数 f 的二阶导数
second_f=diff(first_f);

init_epsilon=1;  % 初始的迭代误差
x_itator=x_init;


% 迭代次数
num=1;  

% 开始迭代
while   init_epsilon>epsilon
   % 将x导数带入到一介导数中
   % symvar(first_f)  是自变量
   df_value=subs(first_f,symvar(first_f),x_itator);
   % 计算二阶导数
    df_2_value=subs(second_f,symvar(second_f),x_itator);

   % 计算
   if df_2_value==0
      df_2_value=double(second_f);
   end

   if  num>nummaxitator
       break;
   end
   delt_x=-(df_value/df_2_value);
   % x2-x1=delt_x
   init_epsilon=abs(delt_x);

   x_itator=x_itator+delt_x;
   num=num+1;
end
num_itorator=num;
x_ok=x_itator;
%  将最优的x 带入的原理的那个方程中 计算出极值
y_ok = subs(f,symvar(f),x_ok);


end

clc;
clear all;
close all;
syms x;
f =  9.*x.^2-sin(x)-1;
[x_optimization,y,num] = Albert_NewTon_Method(f,2,0.000001,300);
x_optimization = double(x_optimization);
y =vpa(y)
x_optimization
num
x = -10:0.01:10;
ft = 9.*x.^2-sin(x)-1;
figure(1)
plot(x,ft);
hold on;
plot(x_optimization,y,'r*')

三、解方程组

假设有两个方程：

 1、先将方程组联立起来

% 将两个方程联立
function  consider_f=ConsiderFunc(x)
 f(1)=x(1)^2-10*x(1)+x(2)^2+8;
 f(2)=x(1)*x(2)^2+x(1)-10*x(2)+8;
 consider_f=[f(1),f(2)];
end

2、计算偏导数，将x 带入就是雅可比矩阵

% 计算偏导  x1'  x2'.......xn'
function  diff_x=Diff_X(x)
% f1   x1  x2  取偏导
% f2   x1  x2  取偏导
%diff_x=[2*x(1)-10,2*x(2) ; x(2)^2 +1,2*x(1)*x(2)-10 ];
diff_x=[2*x(1)-10,2*x(2); x(2)^2+1,2*x(1)*x(2)-10];
end

clear;clc
format;
x0=[0 0];   % 迭代初始值
eps = 0.00001;  % 定位精度要求
for i = 1:100
    f = double(subs(ConsiderFunc(x0),{'x1' 'x2'},{x0(1) x0(2)}));
    df = double(subs(Diff_X(x0),{'x1' 'x2'},{x0(1) x0(2)}));  % 得到雅克比矩阵
    x = x0 - f/df;
    if(abs(x-x0) < eps)
        break;
    end
    x0 = x; % 更新迭代结果
end
disp('定位坐标：');
x
disp('迭代次数：');
i