R 牛顿商及参数传递求解

4. The derivative of a function $f(x)$ at $x$ can be approximated by the Newton's quotient
$$\frac{f(x+h) - f(x-h)}{2 * h},$$
where $h$ is a small number. Write a function to calculate the Newton's quotient for the function $f(x) = \exp(x)$. The function should take two scalar arguments, $x$ and $h$. Use a default value of $h = 1e-6$.
f<-function(x,h=1e-6){
  f<-(exp(x+h)-exp(x-h))/(2*h)
  return(f)
}
f(1)
exp(1)
f(1)-exp(1)
Test your function at the point $x=1$ using the default value of $h$, and compare to the true value of the derivative $f'(1) = e^1$.
5. A very useful feature in R is the ability to pass a function name as an argument.
Here is an example, where 2 is added to the value of a function, for three different functions $\exp(x)$, $\log(x)$, and $\sin(x)$, at selected points $x$.
test <- function(x, f) {
  output <- f(x) + 2
  return(output)
}
test(0, exp)
test(1, log)
test(0, sin)
test(pi/2,sin)
Modify your function from problem 4 so that you pass in the name of the function for which you want to approximate the derivative. Use the same default value for $h$, and approximate the derivative of $sin(x)$ at $x=\pi/4$, of $\log(x)$ at $x=2$, and of $\exp(x)$ at $x=1$.

基于Monster 组和GPT的调写：
函数 f(x) 在 x处的导数可以近似地表示为牛顿商

img

其中 h是一个小的数。编写一个函数来计算函数 f(x) = \exp(x) 的牛顿商。该函数应该接受两个标量参数 x 和 h。使用默认值 h=1e-6。

f <- function(x, h=1e-6) {
  f <- (exp(x+h) - exp(x-h))/(2*h)
  return(f)
}
f(1)
exp(1)
f(1) - exp(1)

用默认值 h 在 x=1 处测试函数，并将其与导数的真实值 f'(1) = e^1 进行比较。

R 中非常有用的一个特性是能够将函数名作为参数传递。以下是一个示例，在选择的点 x 上，对三个不同的函数 \exp(x)，$log(x) 和 \sin(x)的值加上了2。

test <- function(x, f) {
  output <- f(x) + 2
  return(output)
}
test(0, exp)
test(1, log)
test(0, sin)
test(pi/2, sin)

修改问题4中的函数，使其传递要近似导数的函数的名称。使用相同的默认值 $h$，并近似于在 $x=\pi/4 处的 sin(x)，在 x=2 处的 \log(x) 和在 x=1处的 exp(x) 的导数。

f_newton <- function(f, x, h=1e-6) {
  df <- (f(x+h) - f(x-h))/(2*h)
  return(df)
}
f_newton(sin, pi/4)
f_newton(log, 2)
f_newton(exp, 1)

牛顿商是一种用于求解多项式根的方法。假设我们已知一个n次多项式P(x)，其中系数为a0,a1,a2,...,an，且P(x)有一个根为x0。我们希望通过逐步逼近来求解出x0的值。

定义牛顿商：
f(x) = (x - x0)P1(x) + b0
f(x) = (x - x0)((x - x1)P2(x) + b1) + b0
f(x) = (x - x0)((x - x1)((x - x2)P3(x) + b2) + b1) + b0
...
f(x) = (x - x0)((x - x1)((x - x2)...((x - xn-1)Pn(x) + bn-1) + bn-2)... + b0)

其中P1(x), P2(x), ..., Pn(x)为P(x)的每一次差商，也就是：

P1(x) = (P(x) - a0) / (x - x0)
P2(x) = (P1(x) - a1) / (x - x1)
P3(x) = (P2(x) - a2) / (x - x2)
...
Pn(x) = (Pn-1(x) - an-1) / (x - xn-1)

每次求解牛顿商时，我们希望得出下一个根的估计值xk+1，即：

xk+1 = xk - f(xk) / f'(xk)

其中f(x)为上述的牛顿商，f'(x)为其一阶导数，即f'(x) = P1(x) + (x - x0)P2(x) + (x - x0)(x - x1)P3(x) + ... + (x - x0)(x - x1)...(x - xn-2)Pn(x)。

接下来我们使用R语言实现这个算法。首先定义一个函数来计算P(x)的每一次差商：

diff_coef <- function(x, y) {
  if (length(y) == 1) {
    return(y)
  } else {
    return((diff_coef(tail(x, -1), tail(y, -1)) - diff_coef(head(x, -1), head(y, -1))) / (x[length(x)] - x[1]))
  }
}

然后定义一个函数来计算牛顿商：

newton_quot <- function(x, y, epsilon) {
  iter <- 1
  root <- x[1]
  while (iter <= 1000) {
    fx <- y[1]
    dfx <- 1
    for (i in 2:length(x)) {
      dfx <- dfx * (root - x[i-1])
      fx <- fx * (root - x[i-1]) + y[i]
    }
    root_new <- root - fx / dfx
    if (abs(root_new - root) < epsilon) {
      return(root_new)
    }
    root <- root_new
    iter <- iter + 1
  }
  return(NA)
}

最后，我们可以使用上述函数来求解多项式的根：

x <- c(1, 2, 3, 4, 5)
y <- c(-1, 0, 2, 1, -2)
diff_coef(x, y)
newton_quot(x, y, 1e-8)

这里我们以一个5次多项式为例，使用差商函数求解每一次差商，然后使用牛顿商函数求解多项式的根。其中epsilon是我们设置的精度，表示当根的变化量小于这个值时，停止迭代。

该回答引用ChatGPT GPT-4

运行结果

img

将函数名作为参数传递。我们将使用相同的默认值 h，并分别计算在 $x = \pi/4$ 处的 $sin(x)$ 导数，$x = 2$ 处的 $\log(x)$ 导数和 $x = 1$ 处的 $\exp(x)$ 导数。

首先定义了一个名为 newtons_quotient 的函数，该函数接受 x, f 和 h 作为参数，并计算牛顿商。接下来，我们定义了三个点 $x = \pi/4$，$x = 2$ 和 $x = 1$，并分别计算了在这些点上的 $sin(x)$, $\log(x)$ 和 $\exp(x)$ 的导数。最后，我们打印了计算结果。

代码如下


newtons_quotient <- function(x, f, h = 1e-6) {
  quotient <- (f(x + h) - f(x - h)) / (2 * h)
  return(quotient)
}

x1 <- pi / 4
x2 <- 2
x3 <- 1

result_sin <- newtons_quotient(x1, sin)
result_log <- newtons_quotient(x2, log)
result_exp <- newtons_quotient(x3, exp)

print(paste("sin'(pi/4) =", result_sin))
print(paste("log'(2) =", result_log))
print(paste("exp'(1) =", result_exp))