星期日, 二月 05, 2012

用R语言计算微积分

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))

5 条评论:

  1. 我使用density函数得到了两条曲线a,b,想要求mini(a,b)一下的面积如何求解?更为麻烦的是density的长度还不一样

    回复删除
    回复
    1. 可以试下数值积分的做法。

      删除
    2. 谢谢,我之前知道integer函数可以进行数值积分,但是操作的对象是函数表达式。而我使用density得出的是一串有限长度的点,如何在这种情况下使用数值积分?

      删除
    3. 是integrate,前面笔误写错了

      删除
    4. 这方面我不大懂啊,不过有点也可以做吧,先插值,比如样条,生成连续函数,再数值积分。

      删除