在 R 中使用正割法编写一个函数来查找根

Writing a function to find a root using secant method in R

Newton-Raphson Method的标记不正确,我实际上使用的是正割法;但我没有创建新标签的声誉。

我想我离得不远了。我有一个函数,我命名为 ISO() 来计算流体动力学。它是非线性的,因此要找到方程的根,我需要对其进行数值求解。我选择使用正割法,因为我认为它是我想做的最简单的方法。

我认为我遇到问题的地方是重新分配我用来开始迭代的初始 "guesses"。我正在使用 while 循环,我正在计算迭代次数。

sec <- function(x){

  Number_Of_Iterations = 0 # Starts the counter

  x1 = x
  x2 = 3*x

  while(x2 - x1 > 0.0001){ # This is how I'm trying to converge the points.

    Number_Of_Iterations = Number_Of_Iterations + 1

    # The secant function to determine a new x:

    x_new = x2 - ISO(x1)*(x2 - x1)/(ISO(x2) - ISO(x1))

    if(x2 - x_new > x_new - x1){ # These were the rules that I set to reassign the inital chosen values.

      x2 = x_new
      x1 = x1

    }else{

      x1 = x_new
      x2 = x2

    }

  }
  m_dot = x_new
  m_dot

}

当我 运行 这段代码时,我得到的迭代次数是 1,m_dot 的值等于与手动计算第一次迭代时获得的值不同的值。我已经检查过,我的 ISO() 函数 returns 与我手动计算时的值相同,所以底层函数确实有效,我只是认为我无法获得我的 sec() 函数覆盖我所写的内容。

非常感谢所有帮助。

我能猜到发生了什么:初始点处的函数值具有相同的符号,因此割线根 x_new 在区间 [x1, x2] 之外。然后通过重新分配逻辑,您可以在下一次迭代中得到 x2<x1,因此在间隔长度的测试中,差异 x2-x1 为负,因此小于正阈值。因此在没有找到根的情况下迭代后退出。

最有可能的意图是所有的比较都是为了长度,即差异的绝对值。因此更改为

 while( abs(x2 - x1) > 0.0001){ # This is how I'm trying to converge the points.
    ...
    if( abs(x2 - x_new) > abs(x_new - x1) ){ # These were the rules that I set to reassign the inital chosen values.

顺便说一句:根据您的描述,ISO 功能似乎很昂贵。将他们的评估值存储一次并在需要的地方重用这些值,主要是在

x_new = x2 - y2*(x2 - x1)/(y2 - y1)

(请注意,我将错别字从 y1 更正为 y2,您正在求解斜率的等式 (x2-x_new)/(y2-0) = (x2 - x1)/(y2 - y1) 以求割线根。)

因此在循环开始和内部设置 y1=ISO(x1); y2=ISO(x2); 添加到 x 值的重新分配 y 值的计算。

   x2 = x_new; y2 = ISO(x_new);

各自

   x1 = x_new; y1 = ISO(x_new);

这样,对于任何使用的 x 值,ISO 函数只计算一次。

不幸的是,我的问题不仅仅是编码问题;使用正割法时出现问题。我开始迭代的初步估计使函数分析了一个不正确的范围。

LutzL 给出的建议对回答我的问题有很大帮助,abs() 是必要的,他(或她)提出的将函数从割线公式中取出的建议可能会用得更少计算能力。

这是我的工作代码:

sec <- function(x){

  x1 = x
  x2 = x/10 # I have changed the initial values so that x1 is the larger of the two initial estimates.

  while(abs(x2 - x1) > 0.0000001){ # This is how I'm trying to converge the points.

    # The secant function to determine a new x:

    x_new = x2 - square(x2)*(x2 - x1)/(square(x2) - square(x1))

    # I have changed from ISO(), to a simple quadratic formula, I named square(). This secant function works for all of
    # the functions I have since tested it on.

    if( abs(x2 - x_new) > abs(x_new - x1)){

      x2 = x_new

    }else{

      x1 = x_new

    }

  }

  m_dot = x_new
  m_dot

}

我选择编写自己的割线公式而不是使用包中的内置函数的原因是因为这个函数创建起来非常简单,而且我认为使用自己的代码比依赖包更有价值别人写的。

我希望这对刚接触编码的其他人有所帮助,非常感谢那些提供答案的人,尤其是 LutzL,因为如果没有他们的建议,我认为我无法像我一样快地解决这个问题做了。