R語言編程數學分析重讀微積分微分學原理運用

第R語言編程數學分析重讀微積分微分學原理運用目錄1連續(xù)性2求導3數值導數4差商與牛頓插值5方向導數6全微分7法線8偏導數和邊緣檢測基于偏導數的邊緣檢測Roberts算子其他算子

1連續(xù)性

比如下面這個隨機函數

x=seq(0,0.1,0.01)

y=runif(11,0,1)

plot(x,y)

lines(x,y)

x=seq(0,0.01,0.00001)

y=runif(1001,0,1)

plot(x,y)

無論我們把區(qū)間縮小到什么程度，這種亂糟糟的仿佛要鋪滿整個坐標圖的點的樣式并不會發(fā)生變化。

也就是說，f(x)從左邊趨近于0的時候，f(x)在0處是連續(xù)的，而在右側趨近于0的時候，卻并不連續(xù)。此即左連續(xù)和右連續(xù)的概念。

回到切線的問題，如果曲線y=f(x)在x0?點并不連續(xù)，那么這點顯然沒有唯一的一條切線。

有的時候，盡管滿足了連續(xù)性的要求，也不一定存在切線，比如

y=∣x∣

x=seq(-10,10)

y=abs(x)

plot(x,y,type=l)

在x=0的位置，我們找到的切線要么和左邊重合，要么和右邊重合，也就是說這個函數在x=0處存在兩條切線。

同時也就意味著這一點有兩個斜率，兩個導數。所以，如果把導數定義為某種映射，則一個點只能映射為一個值，所以只能定義這點的導數不存在。

3數值導數

根據導數的定義，當函數的定義域不連續(xù)時，其不連續(xù)處顯然是不存在導數的，但圖形可以欺騙我們的眼睛。

x=seq(-1,1,0.1)

y=sin(x)

y1=cos(x)

xEnd=x+0.1

yEnd=y+y1*0.1

plot(x,y)

for(iinseq_along(x)){

+lines(c(x[i],xEnd[i]),c(y[i],yEnd[i]),col=red)

上圖中，圓圈是對y=sinx進行抽樣的結果，可以理解為不連續(xù)函數；紅線則是在每一個分立的點上，sinx在該點的切線。這兩者看上去如此一致，說明連續(xù)函數的導數在抽樣之后仍然具備一定的數學意義。相應地，不連續(xù)的函數，也應該有類似于導數一樣的存在，從而與連續(xù)函數的導數相對應，此即數值導數。如果回顧導數的定義

x=seq(-5,5,0.1)

y=cos(x)

x1=seq(-5,5,0.1)

end=length(x1)

y1=(sin(x1[2:end])-sin(x1[1:end-1]))/0.1

x5=seq(-5,5,0.5)

end=length(x5)

y5=(sin(x5[2:end])-sin(x5[1:end-1]))/0.5

x10=seq(-5,5,1)

end=length(x10)

y10=(sin(x10[2:end])-sin(x10[1:end-1]))/0.5

plot(x,y,type=l,col=red)

lines(x1[2:length(x1)],y1)

lines(x5[2:length(x5)],y5)

lines(x10[2:length(x10)],y10)

如圖所示

由于我們采用的是后向差分，所以三組差商的值整體右移。此外，隨著h的增大，其誤差也越來越明顯。

對一個函數進行反復求導，即可得到高階導數，可以用數學歸納法的方式記為

4差商與牛頓插值

根據這個表達式，可以通過一個簡單的遞歸程序計算數組的差商

#差商算法,x,y為同等長度的數組

diffDiv-function(x,y){

end=length(x)

ind=2:end#索引

return(

if(end==1)y[1]

else(diffDiv(x[ind],y[ind])

-diffDiv(x[ind-1],y[ind-1]))/(x[end]-x[1])

如果要計算階數為k的差商，只需重復調用diffDiv，

kDiffDiv-function(x,y,k){

len-length(x)

if(lenk)return(0)

d-rep(0,len-k)

for(iin1:(len-k))

d[i]-diffDiv(x[i:(i+k)],y[i:(i+k)])

return(d)

據此，繪制出y=x^5這個函數的差商，

plot(x,y)

k=1

lines(x[1:(end-k)],kDiffDiv(x,y,k),col=red)

k=2

lines(x[1:(end-k)],kDiffDiv(x,y,k),col=green)

k=3

lines(x[1:(end-k)],kDiffDiv(x,y,k),col=blue)

k=4

lines(x[1:(end-k)],kDiffDiv(x,y,k),col=purple)

k=5

lines(x[1:(end-k)],kDiffDiv(x,y,k),col=yellow)

k=6

lines(x[1:(end-k)],kDiffDiv(x,y,k),col=gray)

如圖所示

可見差商與微分在階數上有著一致的趨勢。那么，知道了差商之后，可以做點什么呢？

根據差商定義，可得

polyMulti-function(x){

n=length(x)

if(n2)return(c(-x,1))

omega=rep(0,n+1)

omega[1]=-x[1]

omega[2]=1

for(iin2:n){

omega[2:i]=-x[i]*omega[2:i]*+omega[1:(i-1)]

omega[1]=-x[i]*omega[1]

omega[i+1]=1

return(omega)

Newton-function(x,y){

n=length(x)

N=rep(0,n+1)

N[1]=y[1]

for(iin1:n){

omega=polyMulti(x[1:i])

d=kDiffDiv(x[1:i],y[1:i],i-1)

N[1:i]=N[1:i]+d*omega[1:i]

N[i+1]=d*omega[i+1]

return(N)

驗證一下

x=sort(rnorm(10))

y=x^5+2*x^4

N=Newton(x,y)

Y=y*0

for(iin1:length(Y))

for(jin1:length(N))

Y[i]=Y[i]+N[j]*x[i]^(j-1)

plot(x,y)

lines(x,Y)

可見效果還是不錯的。

5方向導數

現繪制出100個隨機點處x方向的偏導數

x=matrix(data=seq(-5,4.95,0.05),nrow=200,ncol=200)

y=t(x)

z=1-(x^2+y^2)

library(rgl)

persp3d(x,y,z,col=red)

N=1000

x0=rnorm(N)

y0=rnorm(N)

z0=1-(x0^2+y0^2)

x1=x0+3

z1=-2*x0*3+z0

for(iin1:N)

lines3d(c(x0[i],x1[i]),c(y0[i],y0[i]),c(z0[i],z1[i]),col=green)

從觀感上來看，綠線的確是沿著xxx方向。但是這個切線顯然不是唯一的，yyy軸方向顯然存在另一條切線。推而廣之，一旦坐標系旋轉，那么曲面上任意一點處的x和y方向均會發(fā)生變化。

那么曲面是否存在一個只和曲面特征有關，而與坐標系無關的參數？

在解決這個問題之前，最好先找到曲面某點沿著任意方向的導數?；仡檟方向偏導數的定義

如果導數的方向發(fā)生旋轉，則可以寫為

如果寫成矢量形式，則定義梯度

沿這些點梯度方向做一些直線，看看效果如何

x=matrix(data=seq(-5,4.95,0.05),nrow=200,ncol=200)

y=t(x)

z=-(x^2+y^2)

library(rgl)

persp3d(x,y,z,col=red)

theta=seq(-pi,pi,0.01)

x0=cos(theta)

y0=sin(theta)

z0=1-(x0^2+y0^2)

x1=x0*0

y1=y0*0

z1=z0+2

for(iin1:N)

lines3d(c(x0[i],x1[i]),c(y0[i],y1[i]),c(z0[i],z1[i]),col=green)

如圖所示，像一頂漂亮的帽子，在某個投影方向看去，和我們熟知的切線如出一轍。

6全微分

做圖如下

theta=seq(-pi,pi,0.1)

x=cos(theta)

y=sin(theta)

plot(x,y,type=l,col=red)

x1=x*0

y1=(x1-x)/(-2*x)*(-2*y)+y

for(iin1:length(theta))

lines(c(x[i],x1[i]),c(y[i],y1[i]),col=green)

可見，梯度方向對應的是圖形的法線方向。對于二維平面內的曲線而言，其法線方向與切線方向垂直。

回顧偏導數的定義

相應地最大方向導數的方向即為梯度的歸一化

現隨機選擇一些點，來繪制一下這四個方向的向量

library(rgl)

N=1500

x-rnorm(N)

y-rnorm(N)

z-1-x^2-y^2

for(iin1:N){

lines3d(c(x[i],3*x[i]),c(y[i],3*y[i]),c(z[i],z[i]+1),col=red)

if(y[i]0.1)

lines3d(c(x[i],x[i]),c(y[i],y[i]-1/y[i]/2),c(z[i],z[i]+1),col=green)

if(x[i]0.1)

lines3d(c(x[i],x[i]-1/x[i]/2),c(y[i],y[i]),c(z[i],z[i]+1),col=green)

lines3d(c(x[i],x[i]*(1-2*z[i])/(2-2*z[i])),c(y[i],y[i]*(1-2*z[i])/(2-2*z[i])),c(z[i],z[i]+1),col=green)

可以看到，綠線幾乎重新編織了一遍原函數，而紅線則刺破了曲面。

8偏導數和邊緣檢測

基于偏導數的邊緣檢測

灰度圖像是天然的z=f(x,y)函數，盡管以一種差分化的形式存在。其中x,y分別代表圖像坐標系中的坐標z可以表示灰度圖像的灰度值。

那么接下來我們可以觀察一下偏導數作用在圖像上是一個什么效果。圖片當然是最經典的lena

library(imager)

img=load.image(lena.jpg)

dim(img)

[1]51251213

gray=grayscale(img)

par(mfrow=c(1,2))

plot(img)

plot(gray)

可見gray顯然為灰度圖像，從其維度就能看得出來，然后將其變?yōu)槎S的數組，接下來就可以進行求導操作了。

dim(gray)

[1]51251211

mat=array(gray,dim=c(512,512))

mat_x=diff(mat,1)

mat_y=t(diff(t(mat),1))

par(mfrow=c(1,2))

image(mat_x)

image(mat_y)

由于圖像坐標系默認是從上向下為yyy軸，從左向右為xxx軸，所以在我們熟知的坐標系中，圖像是上下顛倒的。而且R語言還非常智能（障）地添加了一層偽彩色，這讓我們更加清晰地看出，對圖像進行差分操作，提取出