牛顿迭代法求1/N和A/B的实现方法和matlab代码-优快云博客

本文链接：https://blog.youkuaiyun.com/hy_520520/article/details/146425536

牛顿迭代法

公式推导

假设有一个函数f(x)，要找到f(x)=0得x值。

选取x0作为r的初始近似值，过点(x0,f(x0))做曲线y=f(x)的切线L。

$y=f(x0)+f'(x0)(x-x0)$

则L与x轴交点的横坐标

$x1=x0-\frac{f(x0)}{f'(x0)}$

称x1为r的一次近似值。过点(x1,f(x))做曲线y=f(x)的切线，并求该切线与x轴交点的横坐标

$x2=x1-\frac{f(x1)}{f'(x1)}$

称x2为r的二次近似值。重复上述过程，得r得近似值序列，其中

${x_{n+1}}=x_{n}-\frac{f(x_{n})}{f'(x_{n})}$

称为r得n+1得近似值。称为牛顿迭代公式。

现在想计算任意数据N得倒数1/N。使得x=1/N，可以重新排列成方程N=1/x。所得函数就是

$f(x)=\frac{1}{x}-N$

函数导数就是

$f(x)=-\frac{1}{x^{2}}$

带入牛顿迭代公式

${x_{n+1}}=x_{n}-\frac{(\frac{1}{x_{n}}-N)}{(-\frac{1}{x_{n}^{2}})}$

${x_{n+1}}=x_{n}+x_{n}^{2}(\frac{1}{x_{n}}-N)$

${x_{n+1}}=x_{n}+x_{n}+x_{n}^{2}N$

${x_{n+1}}=x_{n}(2-N*x_{n})$

通过该公式快速逼近x=1/N。

MATLAB实现

任意输入一个数非零数N分解规范化尾数和指数。

$N=m*2^{e},m\epsilon [0.5,1)$

作用是将1/N转换为计算1/m，再利用指数调整结果。

$\frac{1}{N}=(\frac{1}{m})*2^{-e}$

初始值得范围m∈[0.5,1) 所以真实解得范围为1/m∈(1,2]。

如果输入为负数，可以将负数转换为正数计算最后再计算符号位。

sign_m = sign(m);

整理查找表确定初始值，确定查找表得大小，查找表越大，划分得越细，越容易迭代到正常值。设置M=256。将将区间[0.5,1)均匀划分为M个子区间。如果将0.5均匀分256份得话，在第256份是取得得值是0.5+0.5/256*256=1不满足区间大小。所以这里取0.5-255.5。这样就避免了1得存在。

LUT = 1 ./ (0.5 + (0.5/M) * (0.5:1:M-0.5));

通过查表输入当前输入值N的初始值x0。

% 查找初始估计值
index = floor((m - 0.5) / (0.5/M));
index = max(0, min(M-1, index)); % 确保不越界
x = LUT(index + 1);

然后再带入牛顿迭代公式计算。

% 牛顿迭代（3次迭代）
for k = 1:3
    x = x * (2 - m * x);
end

最后再返回结果。完整的matlab代码如下。

clc
clear
%规范化输入将输入N分解为N=m*2^e形式。其中m属于[0.5,1)，e为整数。
%例如N=3，N=0.75*2^2
%分解后查找初始值
% 规范化N到[0.5,1)区间，计算指数exponent
exponent = 0;
N = 6.66666;
sign_m = sign(N);
m=abs(N);
while m >= 1
    m = m / 2;
    exponent = exponent + 1;
end
while m < 0.5
    m = m * 2;
    exponent = exponent - 1;
end

% 生成查找表（示例用256条目）
M = 256;
LUT = 1 ./ (0.5 + (0.5/M) * (0.5:1:M-0.5));
% 查找初始估计值
index = floor((m - 0.5) / (0.5/M));
index = max(0, min(M-1, index)); % 确保不越界
x = LUT(index + 1);

% 牛顿迭代（3次迭代）
for k = 1:3
    x = x * (2 - m * x);
end

% 调整结果并返回
result = sign_m*x * 2^(-exponent);
disp(result);

仿真测试

测试N=3；

测试N=-456；

A/B的计算方法

要计算A/B，最重要的就是计算A*(1/B)再上述的代码上做简单的修改即可。

clc
clear
%规范化输入将输入N分解为N=m*2^e形式。其中m属于[0.5,1)，e为整数。
%例如N=3，N=0.75*2^2
%分解后查找初始值
% 规范化N到[0.5,1)区间，计算指数exponent
exponent = 0;
A =-234;
B =-123;
sign_m = sign(B);
m=abs(B);
while m >= 1
    m = m / 2;
    exponent = exponent + 1;
end
while m < 0.5
    m = m * 2;
    exponent = exponent - 1;
end

% 生成查找表（示例用256条目）
M = 256;
LUT = 1 ./ (0.5 + (0.5/M) * (0.5:1:M-0.5));
% 查找初始估计值
index = floor((m - 0.5) / (0.5/M));
index = max(0, min(M-1, index)); % 确保不越界
x = LUT(index + 1);

% 牛顿迭代（3次迭代）
for k = 1:3
    x = x * (2 - m * x);
end

% 调整结果并返回
format long;
result = A*sign_m*x * 2^(-exponent);
disp('1/N=');
disp(result);

计算

A =234;
B =-123;

结果为