《R的極客理想——進階開發篇 A》一一1.5 R語言的導數計算

本節書摘來華章計算機出版社《r的極客理想——進階開發篇 a》一書中的第1章，第1.5節，作者：張丹 更多章節内容可以通路雲栖社群“異步社群”公衆号檢視。

問題

如何用r語言進行導數計算？

引言

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

不過，r語言為我打開了一扇高數應用的大門，r語言不僅能友善地實作高等數學的計算，還可以很容易地把一篇論文中的高數公式應用于産品的實踐中。因為r語言我重新學習了高數，讓生活中充滿數學，生活會變得更有意思。本節并不是完整的高數計算手冊，僅介紹了導數計算和偏導數計算的r語言實作。

1.5.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=x3的導數，根據導數計算公式，用于手動計算的變形結果為y'=3x2，當x=1時，y'=3，當x=2時，y'=12。

本節的系統環境是：

windows 7 64bit

r: 3.1.1 x86_64-w64-mingw32/x64 (64-bit)

用r語言程式實作導數計算，代碼如下。

dx <- deriv(y ~ x^3, "x") ; dx # 生成導數公式 expression({

}) mode(dx) # 檢視dx變量類型

[1] "expression"

x<-1:2 # 給自變量x指派 eval(dx) # 運作求導計算

[1] 1 8 # 原函數的計算結果

attr(,"gradient") # 使用梯度下降法，導函數的計算結果

[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=π時，y'=-1，當x=4π時，y'=1。

dx <- deriv(y ~ sin(x), "x", func= true) ; dx # 生成導數公式的調用函數 function (x)

{

}

mode(dx) # 檢查dx的類型 [1] "function" dx(c(pi,4*pi)) # 以參數作為自變量，進行函數調用

[1] 1.224606e-16 -4.898425e-16

attr(,"gradient")

[1,] -1 # x=pi,dx=cos(pi)=-1

[2,] 1 # x=4pi,dx=cos(4pi)=1

1.5.2　初等函數的導數公式

對基本的初等函數求導數，通過導數計算公式是可以直接手動完成計算的。下面為一進制初等函數的導數計算公式，其中y是原函數，x是y函數的自變量，y'是y函數的導函數。c,n,a為常數，ln表示以自然常數e為底的對數。

函數 原函數 導函數

常數函數 y=c y'=0

幂函數 y=xn y'=nxn-1

指數函數 y=ax y'=axlna

對數函數 y=logax y'=1/(xln(a)) (a>0,且a!=1,x>0)

正弦函數 y=sin(x) y'=cos(x)

餘弦函數 y=cos(x) y'=-sin(x)

正切函數 y=tan(x) y'=sec2(x)=1/cos2(x)

餘切函數 y=cot(x) y'=-csc2(x)=1/sin2(x)

正割函數 y=sec(x) y'=sec(x)*tan(x)

餘割函數 y=csc(x) y'=-csc(x)*cot(x)

反正弦函數 y=arcsin(x) y'=1/ ()

反餘弦函數 y=arccos(x) y'=-1/ ()

反正切函數 y=arctan(x) y'=1/(1+x2)

反餘切函數 y=arccot(x) y'=-1/(1+x2)

反正割函數 y=arcsec(x) y'=

反餘割函數 y=arccsc(x) y'=

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

1.?常數函數

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

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

[1,] 10 # 導函數計算結果y'=10*1=10

2.?幂函數

計算y=x4函數的導數，根據導數計算公式，用于手動計算的變形結果為y'=4x3，當x=2時，y'=32。

dx<-deriv(y~x^4,"x",func = true) dx(2) [1] 16

[1,] 32 # 導函數計算結果y'=4x^3=42^3=32

3.?指數函數

計算y=4x函數的導數，根據導數計算公式，用于手動計算的變形結果為y'=4xln(4)，當x=2時，y'=22.18071。

0> dx<-deriv(y~4^x ,"x",func = true)

[1,] 22.18071 # 導函數計算結果y'=4^xlog(4)=42^3=22.18071

計算y=ex函數的導數，根據導數計算公式，用于手動計算的變形結果為y'=ex，當x=2時，y'=y=7.389?056。

dx<-deriv(y~exp(1)^x ,"x",func = true) [1] 7.389056

[1,] 7.389056 # 導函數計算結果y'=exp(1)^x=exp(1)^2=7.389056

4.?對數函數

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

dx<-deriv(y~log(x),"x",func = true) [1] 0.6931472

[1,] 0.5 # 導函數計算結果y'=1/x=1/2=0.5

計算y=log2x函數的導數，根據導數計算公式，用于手動計算的變形結果為y'=

1/xlna，當x=3時，y'=0.480?898?3。但用r語言程式設計時，隻能計算以自然常數為底的對數的導數，對于原函數不是以自然常數為底的對數，首先要變換成以自然常數為底的對數再進行導數計算，根據對數的換底公式，把以2為底的對數轉換為以自然常數為底的對數y=log2x=，

dx<-deriv(y~log(x)/log(2),"x",func = true) dx(3) [1] 1.584963

[1,] 0.4808983 # 導函數計算結果y'=1/(xlog(2)=1/(3log(2)=0.4808983

5.?正弦函數

計算y=sin(x)函數的導數，根據導數計算公式，用于手動計算的變形結果為y'=cos(x)，當x=π時，y'=-1。

dx<-deriv(y~sin(x),"x",func = true) dx(pi) [1] 1.224606e-16

[1,] -1 # 導函數計算結果y'=cos(x)=cos(pi)=-1

6.?餘弦函數

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

dx<-deriv(y~cos(x),"x",func = true) dx(pi/2) [1] 6.123032e-17

[1,] -1 # 導函數計算結果y'=-sin(x)=-sin(pi/2)=-1

7.?正切函數

計算y=tan(x)函數的導數，根據導數計算公式，用于手動計算的變形結果為y'=sec2(x)= 1/cos2(x)，當x=π/6時，y'=1.333?333。

dx<-deriv(y~tan(x),"x",func = true) dx(pi/6) [1] 0.5773503

[1,] 1.333333 # 導函數計算結果y'=1/cos(x)^2=1/cos(pi/6)^2=1.333333

8.?餘切函數

計算y=cot(x)函數的導數，由于r語言沒有cot()函數，是以根據三角公式我們手動變形原函數為y=cot(x)=1/tan(x)後再進行導數計算，根據導數計算公式，用于手動計算的變形結果為y'=-csc2(x)=-1/sin2(x)，當x=π/6時，y'=-4。

dx<-deriv(y~1/tan(x),"x",func = true) [1] 1.732051

[1,] -4 # 導函數計算結果y'=-1/sin(x)^2=-1/sin(pi/6)^2=-4

9.反正弦函數

計算y=arcsin(x)函數的導數，根據導數計算公式，用于手動計算的變形結果為y'=1/，當x=π/6時，y'=1.173?757。

dx<-deriv(y~asin(x),"x",func = true) [1] 0.5510696

[1,] 1.173757 # 導函數計算結果y'=1/sqrt(1-x^2)=1/sqrt(1-(pi/6)^2)=1.173757

10.?反餘弦函數

計算y=arccos(x)函數的導數，根據導數計算公式，用于手動計算的變形結果為y'=

-1/，當x=π/8時，y'=-1.087?35。

dx<-deriv(y~acos(x),"x",func = true) dx(pi/8) [1] 1.167232

[1,] -1.08735 # 導函數計算結果y'=-1/sqrt(1-x^2)=-1/sqrt(1-(pi/8)^2)=-1.08735

11.?反正切函數

計算y=arctan(x)函數的導數，根據導數計算公式，用于手動計算的變形結果為y'=

1/(1+x2)，當x=π/6時，y'=0.784?833?5。

dx<-deriv(y~atan(x),"x",func = true) [1] 0.4823479

[1,] 0.7848335 # 導函數計算結果y'= 1/(1+x^2) = 1/(1+(pi/6)^2)=0.7848335

1.5.3　二階導數計算

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

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

下面計算函數y=sin(ax)的二階導數y''，其中a為常數。根據導數計算公式，用手動計算的變形結果，一階導數為y'=acos(ax)，對y'再求導公式變形為，y''=-a2sin(ax)

用r語言進行程式實作

a<-2 # 設定a的值 dx<-deriv(y~sin(a*x),"x",func = true) # 生成一階導數公式 dx(pi/3) # 計算一階導數 [1] 0.8660254

[1,] -1 # 導函數計算結果y'= acos(ax)=2cos(2pi/3)=-1

dx<-deriv(y~acos(ax),"x",func = true) # 對一階導函數求導 dx(pi/3) [1] -1

[1,] -3.464102 # 導函數計算結果y'= -a^2sin(ax)=-2^2sin(2pi/3)=-3.464102

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

dx<-deriv3(y~sin(a*x),"x",func = true) # 生成二階導數公式 dx(pi/3) # 計算導數

[1,] -1 # 一階導數結果

attr(,"hessian")

, , x

[1,] -3.464102 # 二階導數結果

我們再計算另外一個二階導數，計算y=ax4+bx3+x2+x+c，其中a, b, c為常數，a=2, b=1, c=3，根據導數計算公式，手動計算的變形結果，一階導數為y'=(2x4+x3+x2+x+3'?)=?8x3+3x2+

2x+1，當x=2時，y'=81，對y'再求導公式變形為，y''=24x2+6x+2，當x=2時，y''=110。

dx<-deriv3(y~ax^4+bx^3+x^2+x+c,"x",func=function(x,a=2,b=1,c=3){})

dx(2) [1] 49

[1,] 81 # 一階導數結果

[1,] 110 # 二階導數結果

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

1.5.4　偏導數計算

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

在xoy平面内，當動點由p(x0, y0)沿不同方向變化時，函數f(x, y)的變化快慢一般來說是不同的，是以就需要研究f(x, y)在點(x0, y0)處沿不同方向的變化率。在這裡我們隻學習函數f(x, y)在xoy平面沿着平行于x軸和平行于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)=2x2+y+3xy2的偏導數，由于二進制函數曲面上每一點都有無窮多條切線，描述這個函數的導數就會相當困難。如果讓其中的一個變量y取值為常數，那麼就可以求出關于另一個自變量x的偏導數了，即?f/?x。

下面我們分别對x, y兩個自變量求偏導數，設變量y為常數，計算x的偏導數?f/?x= 4x+3y2，當x=1, y=1時，x的偏導數?f/?x=4x+3y2=7。設變量x為常數，計算y的偏導數?f/?y= 1+6xy，當x=1, y=1時，y的偏導數?f/?x=1+6xy=7。r語言程式實作如下。

fxy = expression(2x^2+y+3x*y^2) # 二進制函數公式 dxy = deriv(fxy, c("x", "y"), func = true) dxy function (x, y)

dxy(1,1) # 設定自變量 [1] 6

[1,] 7 7

偏導數的程式計算結果與手動計算結果是一緻的。

下面我們再求一個複雜函數的偏導數，計算一個二進制函數f(x, y)=x?y+exy+x2-2xy+ y3+sin(xy)在點(1, 3)和點(0, 0)的偏導數。r語言程式實作如下。

fxy = expression(x^y + exp(x y) + x^2 - 2 x y + y^3 + sin(xy)) dxy(1,3) # 設定自變量 [1] 43.22666

[1,] 56.28663 44.09554 # 計算結果，x的偏導數為56.28663，y的偏導數為 44.09554

dxy(0,0) [1] 2

[1,] nan -inf # 計算結果，x的偏導數無意義，y的偏導數負無窮大

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

本節我們學習了用r語言做高等數學的導數計算，真的是非常友善，這下更有動力學習高數了。