R語言的導數計算（轉）

轉自：http://blog.fens.me/r-math-derivative/

前言

高等數學是每個大學生都要學習的一門數學基礎課，同時也可能是考完試后最容易忘記的一門知識。
我在學習高數的時候絞盡腦汁，但始終都不知道為何而學。生活和工作基本用不到，就算是在計算機行業和金融行業，能直接用到高數的地方也少之又少，學術和實際應用真是相差太遠了。

不過，R語言為我打開了一道高數應用的大門，R語言不僅能方便地實現高等數學的計算，還可以很容易地把一篇論文中的高數公式應用於產品的實踐中。
因為R語言我重新學習了高數，讓生活中充滿數學，生活會變得更有意思。

本節並不是完整的高數計算手冊，僅介紹了導數計算和偏導數計算的R語言實現。

目錄

1. 導數計算
2. 初等函數的導數公式
3. 二階導數計算
4. 偏導數計算

1. 導數計算

導數(Derivative)是微分學的基本概念，用於計算函數的極值。導數的定義為，當函數y=f(x)在x0的某個領域內有定義，當自變量x在x0處取得增加Δx(點x0+Δx仍在該鄰域內)時，相應的函數取得增量Δy=f(x0+Δx)-f(x0)；如果Δy與Δx之比當Δx趨於0時的極限存在，則稱函數y=f(x)在點x0處可導，並稱這個極限為函數y=f(x)在點x0處的導數，記為f`(x0)，即

也記作 y’|x=x0 ，dy/dx|x=x0 或 df(x)/dx|x=x0。

通過R語言可以使用deriv()函數直接進行導數的計算，比如要計算 y=x^3 的導數，根據導數計算公式，用於手動計算的變形結果為 y’=3x^2，當x=1時，y’=3,當x=2時，y’=12。

本節的系統環境

• Win7 64bit
• R: 3.1.1 x86_64-w64-mingw32/x64 (64-bit)

用R語言程序實現，代碼如下。

> dx <- deriv(y ~ x^3, "x") ; dx  # 生成導數公式
expression({
    .value <- x^3
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .grad[, "x"] <- 3 * x^2
    attr(.value, "gradient") <- .grad
    .value
})
> mode(dx)                        # 查看dx變量類型
[1] "expression"
> x<-1:2                          # 給自變量x賦值
> eval(dx)                        # 運行求導計算
[1] 1 8                           # 原函數的計算結果
attr(,"gradient")                 # 使用梯度下降法，導函數的計算結果
      x
[1,]  3                           # x=1,dx=3*1^2=3
[2,] 12                           # x=2,dx=3*2^2=12

用R語言程序計算的結果，與我們手動計算的結果是一致的。但計算過程其實是有很大區別的，我們手動計算時是通過給定的導數計算公式，變成后完成的計算。而用計算機程序計算時，是使用梯度下降法來計算一階導數，是一種最優化的近似算法。對於手動計算導數時，如果函數比較復雜而且比較難應用可變形的公式，那么手動計算就會有非常大的困難，而計算機程序的方法是一般地導數計算方法，不會受到公式難於變形的影響。

我們使用deriv(expr, name)函數時通常要傳2個參數，第一參數expr就是原函數公式，用~號來分隔公式的兩邊，第二參數name用於指定函數的自變量。deriv()函數會返回一個表達式expression類型變量，再用eval()函數運行這個表達式得到就可得到計算結果，如上面的代碼實現。

如果希望以函數的形式調用計算公式，那么你還需要傳第三個參數func，並讓func參數為TRUE，參考下面的代碼實現。

計算正弦函數y=sin(x)的導數，根據導數計算公式，用於手動計算的變形結果為 y’=cos(x)，當x=pi時，y’=-1，當x=4*pi時，y’=1，其中pi=π表示圓周率。

> dx <- deriv(y ~ sin(x), "x", func= TRUE) ; dx      # 生成導數公式的調用函數
function (x)
{
    .value <- sin(x)
    .grad <- array(0, c(length(.value), 1L), list(NULL, c("x")))
    .grad[, "x"] <- cos(x)
    attr(.value, "gradient") <- .grad
    .value
}
> mode(dx)                               # 檢查dx的類型
[1] "function"
> dx(c(pi,4*pi))                         # 以參數作為自變量，進行函數調用
[1]  1.224606e-16 -4.898425e-16
attr(,"gradient")
      x                                  # 導函數的計算結果
[1,] -1                                  # x=pi,dx=cos(pi)=-1
[2,]  1                                  # x=4*pi,dx=cos(4*pi)=1

2. 初等函數的導數公式

對於基本的初等函數求導數，通過導數計算公式是可以直接手動完成計算的，下面為一元初等函數的導數計算公式。

函數          原函數            導函數
常數函數      y=C               y'=0
冪函數        y=x^n             y'=n*x^(n-1)
指數函數      y=a^x             y'=a^x*ln(a)
              y=exp(1)^x        y'=exp(1)^x
對數函數      y=log(x,base=a)   y'=1/(x*ln(a)) (a>0，且a!=1,x>0)
              y=ln(x)           y'=1/x
正弦函數      y=sin(x)          y'=cos(x)
余弦函數      y=cos(x)          y'=-sin(x)
正切函數      y=tan(x)          y'=sec(x)^2=1/cos(x)^2
余切函數      y=cot(x)          y'=-csc(x)^2=1/sin(x)^2
正割函數      y=sec(x)          y'=sec(x)*tan(x)
余割函數      y=csc(x)          y'=-csc(x)*cot(x)
反正弦函數    y=arcsin(x)       y'=1/sqrt(1-x^2)
反余弦函數    y=arccos(x)       y'=-1/sqrt(1-x^2)
反正切函數    y=arctan(x)       y'=1/(1+x^2)
反余切函數    y=arccot(x)       y'=-1/(1+x^2)
反正割函數    y=arcsec(x)       y'=1/abs(x)*(x^2-1)
反余割函數    y=arccsc(x)       y'=-1/abs(x)*(x^2-1)

公式的注釋：

• y是原函數，x是y函數的自變量，y’是y函數的導函數。
• C,n,a為常數。
• ln表示以自然常數e為底的對數
• exp(1)表示自然常數e
• log(x,base=a)表示，以常數a為底的對數
• sqrt表示開平方
• abs表示絕對值
• 正割函數sec，計算方法為 sec=1/cos(x)
• 余割函數csc，計算方法為 csc=1/sin(x)
• 余切函數cot，計算方法為 cot=1/tan(x)

注： 以上公式不完全匹配於R語言函數

接下來，我們分別對這些一元初等函數進行一階導數的計算。設y為原函數，x是y函數的自變量，且只有一個自變量。

常數函數

計算 y=3+10*x 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=0+10*x ，常數項3的導數為0，當x=1時，y’=10。

> dx<-deriv(y~ 3+10*x,"x",func = TRUE)  # 以函數形式生成導數公式
> dx(1)                                 # 傳入自變量，並計算
[1] 13                                  # 原函數計算結果y=3+10*1=13
attr(,"gradient")
      x
[1,] 10         # 導函數計算結果y'=10*1=10

冪函數

計算 y=x^4 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=4*x^3，當x=2時，y’=32。

> dx<-deriv(y~x^4,"x",func = TRUE)
> dx(2)
[1] 16
attr(,"gradient")
      x
[1,] 32         # 導函數計算結果y'=4*x^3=4*2^3=32

指數函數

計算 y=4^x 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=4^x*ln(4)，當x=2時，y’=22.18071。

> dx<-deriv(y~4^x ,"x",func = TRUE)
> dx(2)
[1] 16
attr(,"gradient")
            x
[1,] 22.18071   # 導函數計算結果y'=4^x*log(4)=4*2^3=22.18071

計算 y=exp(1)^x 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=exp(1)^x，當x=2時，y’=y=7.389056。

> dx<-deriv(y~exp(1)^x ,"x",func = TRUE)
> dx(2)
[1] 7.389056
attr(,"gradient")
            x
[1,] 7.389056   # 導函數計算結果y'=exp(1)^x=exp(1)^2=7.389056

對數函數

計算 y=ln(x) 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=1/x，當x=2時，y’=0.5。

> dx<-deriv(y~log(x),"x",func = TRUE)
> dx(2)
[1] 0.6931472
attr(,"gradient")
       x
[1,] 0.5   # 導函數計算結果y'=1/x=1/2=0.5

計算 y=log2(x) 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=1/(x*log(2))，當x=3時，y’=0.4808983。

但用R語言編程時，只能計算以自然常數為底的對數的導數，對於原函數不是以自然常數為底的對數，首先要變換成以自然常數為底的對數再進行導數計算，根據對數的換底公式，把以2為底的對數轉換為以自然常數為底的對數 y=log2(x)=log(x)/log(2)，

> dx<-deriv(y~log(x)/log(2),"x",func = TRUE)
> dx(3)
[1] 1.584963
attr(,"gradient")
             x
[1,] 0.4808983         # 導函數計算結果y'=1/(x*log(2)=1/(3*log(2)=0.4808983

正弦函數

計算 y=sin(x) 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=cos(x)，當x=pi時，y’=-1，其中pi=π表示圓周率。

> dx<-deriv(y~sin(x),"x",func = TRUE)
> dx(pi)
[1] 1.224606e-16
attr(,"gradient")
      x
[1,] -1   # 導函數計算結果y'=cos(x)=cos(pi)=-1

余弦函數

計算 y=cos(x) 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=-sin(x)，當x=pi/2時，y’=-1。

> dx<-deriv(y~cos(x),"x",func = TRUE)
> dx(pi/2)
[1] 6.123032e-17
attr(,"gradient")
      x
[1,] -1  # 導函數計算結果y'=-sin(x)=-sin(pi/2)=-1

正切函數

計算 y=tan(x) 函數的導數，根據導數計算公式，用於手動計算的變形結果為 y’=sec(x)^2=1/cos(x)^2，當x=pi/6時，y’=1.333333。

> dx<-deriv(y~tan(x),"x",func = TRUE)
> dx(pi/6)
[1] 0.5773503
attr(,"gradient")
            x
[1,] 1.333333  # 導函數計算結果y'=1/cos(x)^2=1/cos(pi/6)^2=1.333333

余切函數

計算 y=cot(x) 函數的導數，由於R語言沒有cot()函數，所以根據三角公式我們動手變形原函數為y=cot(x)=1/tan(x)后再進行導數計算，根據導數計算公式，用於手動計算的變形結果為y’=-csc(x)^2=-1/sin(x)^2，當x=pi/6時，y’=-4。

> dx<-deriv(y~1/tan(x),"x",func = TRUE)
> dx(pi/6)
[1] 1.732051
attr(,"gradient")
      x
[1,] -4  # 導函數計算結果y'=-1/sin(x)^2=-1/sin(pi/6)^2=-4

反正弦函數

計算 y=asin(x) 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=1/sqrt(1-x^2)，當x=pi/6時，y’=1.173757。

> dx<-deriv(y~asin(x),"x",func = TRUE)
> dx(pi/6)
[1] 0.5510696
attr(,"gradient")
            x
[1,] 1.173757  # 導函數計算結果y'=1/sqrt(1-x^2)=1/sqrt(1-(pi/6)^2)=1.173757

反余弦函數

計算 y=acos(x) 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=-1/sqrt(1-x^2)，當x=pi/8時，y’=-1.08735。

> dx<-deriv(y~acos(x),"x",func = TRUE)
> dx(pi/8)
[1] 1.167232
attr(,"gradient")
            x
[1,] -1.08735    # 導函數計算結果y'=-1/sqrt(1-x^2)=-1/sqrt(1-(pi/8)^2)=-1.08735

反正切函數

計算 y=atan(x) 函數的導數，根據導數計算公式，用於手動計算的變形結果為y’=1/(1+x^2)，當x=pi/6時，y’=0.7848335。

> dx<-deriv(y~atan(x),"x",func = TRUE)
> dx(pi/6)
[1] 0.4823479
attr(,"gradient")
             x
[1,] 0.7848335   # 導函數計算結果y'= 1/(1+x^2) = 1/(1+(pi/6)^2)=0.7848335

3. 二階導數計算

當我們對一個函數進行多次接連的求導計算，會形成高階導數。
一般的，函數y=f(x)的導數y’=f’(x)仍然是x的函數，我們就把y’=f’(x)的導數叫做函數y=f(x)的二階導數，記作y”，即

一階導數的導數叫做二階導數，二階導數的導數叫做三階導數，N-1階導數的導數叫做N階導數，習慣上把二階以上的導數稱之為高階導數，

比如，計算 y=sin(a*x) 函數的二階導數導數y”，其中a為常數，根據導數計算公式，用於手動計算的變形結果為一階導數為y’=a*cos(a*x)，對y’再求導公式變形為，y”=-a^2*sin(a*x)

用R語言進行程序實現

> a<-2                 # 設置a的值
> dx<-deriv(y~sin(a*x),"x",func = TRUE) # 生成一階導數公式
> dx(pi/3)                              # 計算一階導數
[1] 0.8660254
attr(,"gradient")
      x
[1,] -1     # 導函數計算結果y'= a*cos(a*x)=2*cos(2*pi/3)=-1

> dx<-deriv(y~a*cos(a*x),"x",func = TRUE)    # 對一階導函數求導
> dx(pi/3)
[1] -1
attr(,"gradient")
             x
[1,] -3.464102     # 導函數計算結果y'= -a^2*sin(a*x)=-2^2*sin(2*pi/3)=-3.464102

上面二階導數的計算，我們是動手划分為兩次求導進行計算的，利用deriv3()函數其實合並成一步計算。

> dx<-deriv3(y~sin(a*x),"x",func = TRUE)  # 生成二階導數公式
> dx(pi/3)                                # 計算導數
[1] 0.8660254
attr(,"gradient")
      x
[1,] -1         # 一階導數結果
attr(,"hessian")
, , x
             x
[1,] -3.464102  # 二階導數結果

我們再計算另外一個二階導數，計算y=a*x^4+b*x^3+x^2+x+c，其中a,b,c為常數a=2,b=1,c=3，
根據導數計算公式，用於手動計算的變形結果為一階導數為y’=2*x^4+x^3+x^2+x+3=4*2*x^3+3*x^2+2*x+1,當x=2時，y’=81，
對y’再求導公式變形為，y”=3*4*2*x^2+2*3*x+2，當x=2時，y”=110。

> dx<-deriv3(y~a*x^4+b*x^3+x^2+x+c,"x",func=function(x,a=2,b=1,c=3){})  # 通過func參數，指定常數值
> dx(2)
[1] 49
attr(,"gradient")
      x
[1,] 81           # 一階導數結果
attr(,"hessian")
, , x
       x
[1,] 110          # 二階導數結果

這樣就直接完成了二階導數的計算，在R語言中二階導數是可以直接求出的，想計算更高階的導數就需要其他的數學工具包了。

4. 偏導數計算

在一元函數中，我們已經知道導數就是函數的變化率。對於二元函數我們同樣要研究它的“變化率”。然而，由於自變量多了一個，情況就要復雜的多。在數學中，一個多變量的函數的偏導數，就是它關於其中一個變量的導數而保持其他變量恆定（相對於全導數，在其中所有變量都允許變化）。
偏導數的算子符號為:∂。記作∂f/∂x 或者 f’x。偏導數反映的是函數沿坐標軸正方向的變化率，在向量分析和微分幾何中是很有用的。

在xOy平面內，當動點由P(x0,y0)沿不同方向變化時，函數f(x,y)的變化快慢一般說來是不同的，因此就需要研究f(x,y)在(x0,y0)點處沿不同方向的變化率。在這里我們只學習函數f(x,y)在x0y平面沿着平行於x0y軸和平行於y軸兩個特殊方位變動時，f(x,y)的變化率。

x方向的偏導：
設有二元函數z=f(x,y)，點(x0,y0)是其定義域D內一點.把y固定在y0而讓x在x0有增量△x，相應地函數z=f(x,y)有增量(稱為對x的偏增量)△z=f(x0+△x,y0)-f(x0,y0)。如果△z與△x之比當△x→0時的極限存在，那么此極限值稱為函數z=f(x,y)在(x0,y0)處對x的偏導數(partial derivative)。記作f’x(x0,y0)。

y方向的偏導：
函數z=f(x,y)在(x0,y0)處對x的偏導數，實際上就是把y固定在y0看成常數后，一元函數z=f(x,y0)在x0處的導數。同樣，把x固定在x0,讓y有增量△y,如果極限存在那么此極限稱為函數z=(x,y)在(x0,y0)處對y的偏導數。記作f’y(x0,y0)

同樣地，我們可以通過R語言的 deriv()函數進行偏導數的計算。下面我們計算一個二元函數f(x,y)=2*x^2+y+3*x*y^2的偏導數，由於二元函數曲面上每一點都有無窮多條切線，描述這個函數的導數就會相當困難。如果讓其中的一個變量y取值為常數，那么就可以求出關於另一個自變量x的偏導數了，即∂f/∂x。

下面我們分別對x,y兩個自變量求偏導數，設變量y為常數，計算x的偏導數∂f/∂x=4*x+3*y^2，當x=1,y=1時，x的偏導數∂f/∂x=4*x+3*y^2=7。設變量x為常數，計算y的偏導數∂f/∂y=1+6*x*y，當x=1,y=1時，y的偏導數∂f/∂x=1+6*x*y=7。

用R語言程序實現。

> fxy = expression(2*x^2+y+3*x*y^2)     # 二元函數公式
> dxy = deriv(fxy, c("x", "y"), func = TRUE)
> dxy
function (x, y)
{
    .expr4 <- 3 * x
    .expr5 <- y^2
    .value <- 2 * x^2 + y + .expr4 * .expr5
    .grad <- array(0, c(length(.value), 2L), list(NULL, c("x","y")))
    .grad[, "x"] <- 2 * (2 * x) + 3 * .expr5
    .grad[, "y"] <- 1 + .expr4 * (2 * y)
    attr(.value, "gradient") <- .grad
    .value
}
> dxy(1,1)                          # 設置自變量
[1] 6
attr(,"gradient")
     x y                            # 計算結果，x的偏導數為7，y的偏導數為7
[1,] 7 7

偏導數的程序計算結果與手動計算結果是一致的。下面我們再求一個復雜函數偏導數，計算一個二元函數f(x,y)=x^y + exp(x * y) + x^2 – 2 * x * y + y^3 + sin(x*y)在點(1,3)和點(0,0)的偏導數。

R語言程序實現。

> fxy = expression(x^y + exp(x * y) + x^2 - 2 * x * y + y^3 + sin(x*y))
> dxy = deriv(fxy, c("x", "y"), func = TRUE)
> dxy(1,3)        # 設置自變量
[1] 43.22666
attr(,"gradient")
            x        y
[1,] 56.28663 44.09554        # 計算結果，x的偏導數為56.28663，y的偏導數為 44.09554
> dxy(0,0)
[1] 2
attr(,"gradient")
       x    y
[1,] NaN -Inf                 # 計算結果，x的偏導數無意義，y的偏導數負無窮大

對於計算的結果，有異議的同學，可以嘗試動手計算。

本文我們掌握了R語言對於高等數學的導數計算方法，真的是非常方便，這下更有動力學習高數了。