Matlab中的Newton-Raphson方法

Question

我是matlab的新手,我需要创建一个函数,该函数使用起始近似x = a进行N次迭代的Newton-Raphson方法.这种起始近似不算作交互,另一个要求是需要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技术的定义:

blah http://web.mit.edu/10.001/Web/Course_Notes/NLAE/equation6.gif

对于下一次迭代,您使用上一次迭代的值.你正在做的是使用循环计数器并将其替换为你的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的这篇文章,你会发现使用这些变量作为循环索引实际上更慢:在Matlab中使用i和j作为变量

在任何情况下,为了测试这个,假设我们的函数是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.