用R语言计算微积分

This post was kindly contributed by 数据科学与R语言 - go there to comment and to read the full post.

R语言是一个优秀的统计计算环境,但它能计算微积分吗?回答是:这个可以有。虽然在这方面R比不上Mathmatica或是Matlab,不太复杂的问题仍然能解决。以下图为例,我们希望计算两个简单问题,一个是sin函数在x=4处的导数。另一个问题是计算sin(x)在0到pi之间的曲线下面积。
第一个问题用D运算子得到符号运算结果,然后代入数值即可得到最后结果
dx <- D(expression(sin(x)),'x')
x <- 4
slope <- eval(dx)

第二个问题可以用integrate函数直接求出数值积分
integrate(function(x) sin(x),0,pi)
数值积分的另一种算法是蒙特卡罗仿真,得到的结果相当接近
x <- runif(100000,min=0,max=pi)
y <- runif(100000)
pi*sum(y<sin(x))/100000
上面的图形是用ggplot2扩展包绘制的,用到的代码如下:
library(ggplot2)
intercept <- sin(4)-slope*4
x <- seq(from=0,to=2*pi,by=0.01)
y <- sin(x)
p <- ggplot(data.frame(x,y),aes(x,y))
p + geom_area(fill=alpha('blue',0.3))+
geom_abline(intercept=intercept,slope=slope,linetype=2)+
scale_x_continuous(breaks=c(0,pi,2*pi),labels=c('0',expression(pi),expression(2*pi)))+
geom_text(parse=T,aes(x=pi/2,y=0.3,label='integral(sin(x)*dx, 0, pi)'))+
geom_line()+
geom_point(aes(x=4,y=sin(4)),size=5,colour=alpha('red',0.5))