matlab - Matlab 中的 Newton-Raphson 方法

Question

我是 matlab 的新手，我需要创建一个函数，以起始近似值 x = a 执行 Newton-Raphson 方法的 n 次迭代。此起始近似值不算作交互，另一个要求是需要 for 循环。我查看了发布的其他类似问题，但就我而言，我不想使用 while 循环。

我的输入应该是这样的:

mynewton(f,a,n) which takes three inputs: 
f: A function handle for a function of x.
a: A real number.
n: A positive integer.

到目前为止，这是我的代码。

function r=mynewton(f,a,n)
syms x;
z=f(x);
y=a;
for i=1:n    
    y(i+1)=y(i)-(z(i)/diff(z(i)));
end
r=y
end

当我尝试调用该函数时，我收到一条错误消息:

Error in MuPAD command: DOUBLE cannot convert the input expression into a double    array.
If the input expression contains a symbolic variable, use the VPA function instead.
Error in mynewton (line 6)
y(i+1)=y(i)-(z(i)/diff(z(i)));

问题是我该如何使用这个VPA功能？当然，我的代码的其余部分可能也不是 100% 正确，但是如果能解决 vpa 问题或修复我的代码的其他部分，我们将不胜感激。

谢谢!

Answer 1

您的 Newton-Raphson 技术有两点不太正确……但肯定可以解决!在我们修复此问题后，您所说的 VPA 错误就不再存在了。

---

错误 #1 - 迭代更新

第一个是迭代本身。记忆一下 Newton-Raphson 技术的定义:

_{(来源:mit.edu)}

对于下一次迭代，您使用上一次迭代的值。您所做的是使用循环计数器并将其代入您的f(x)，这是不正确的。它必须是前一次迭代的值。

错误 #2 - 将符号值与数值混合

如果您看一下您是如何编写函数代码的，您会以符号方式定义您的函数，但您正在尝试将数字 值替换到您的函数中。不幸的是，这不适用于 MATLAB。如果你真的想替换值，你必须使用 subs .这将为您替换一个实际值作为 x 的函数或您的函数使用的任何自变量。执行此操作后，您的值仍然是 sym 类型。您需要将其转换为 double 才能以数字方式使用它。

---

同样为了效率，没必要把y做成数组。只需将其设为在每次迭代时更新自身的单个值即可。综上所述，您的代码已更新为如下所示。请注意，我在循环之前采用了函数的导数来减少您需要进行的计算量。我还拆分了 Newton-Raphson 迭代的分子和分母项，以使事情变得清晰，并使 subs 更容易接受。事不宜迟:

function r = mynewton(f,a,n)
syms x;
z = f(x);
diffZ = diff(z); %// Edit - Include derivative
y = a; %// Initial root

for idx = 1 : n    
    numZ = subs(z,x,y); %// Numerator - Substitute f(x) for f(y)
    denZ = subs(diffZ,x,y); %// Denominator - Substitute for f'(x) for f'(y)
    y = y - double(numZ)/double(denZ); %// Update - Cast to double to get the numerical value
end
r = y; %// Send to output
end

---

请注意，我在循环中将 i 替换为 idx。原因是因为实际上不建议使用 i 或 j 作为循环索引，因为这些字母被保留来表示复数。如果你看看 Shai 发表的这篇文章，你会发现使用这些变量作为循环索引实际上更慢:Using i and j as variables in Matlab

无论如何，为了测试这一点，假设我们的函数是 y = sin(x)，我的初始根是 x0 = 2，有 5 次迭代，我们这样做:

f = @(x) sin(x);
r = mynewton(f, 2, 5)

r =

3.1416

这与我们对 sin(x) 的了解一致，因为 sin(x) 的截距位于 pi 的整数倍处. x0 = 2 位于 pi 附近，因此这确实按我们预期的那样工作。

---

给你的小红包

您的原始代码在每次迭代时将根的值存储在 y 中。如果您真的想这样做，则必须修改您的代码，使其看起来像这样。请记住，我预先分配了 y 以提高效率:

function r = mynewton(f,a,n)
syms x;
z = f(x);
diffZ = diff(z);
y = zeros(1,n+1); %// Pre-allocate output array 
y(1) = a; %// First entry is the initial root

for idx = 1 : n    
    numZ = subs(z,x,y(idx)); %// Remember to use PREVIOUS guess for next guess
    denZ = subs(diffZ,x,y(idx));
    y(idx+1) = y(idx) - double(numZ)/double(denZ); %// Place next guess in right spot  
end
r = y; %// Send to output
end

通过使用与上面完全相同的参数运行此代码，我们得到:

f = @(x) sin(x);
r = mynewton(f, 2, 5)

r =

    2.0000    4.1850    2.4679    3.2662    3.1409    3.1416

r 中的每个值都会告诉您在该特定迭代中对根的猜测。数组的第一个元素是初始猜测(当然)。下一个值是 Newton-Raphson 根每次迭代时的猜测值。请注意，数组的最后一个元素是我们的最终迭代，大致等于 pi。