如何用R求定积分阿

Python025

如何用R求定积分阿,第1张

用R求定积分应该是指积分区间是R吧！我不知道怕我的理解对吗？

如果区间是R的话，就是那些反常积分中的无穷限积分，一般可以用牛顿-莱布尼茨公式，首先把积分函数写成其原函数的形式，然后在根据求极限的方法可以求出在正无穷或负无穷的值就OK啦！

还可以把无穷限积分的积分区间分割成几个小部分在再计算

一般定积分是很灵活的。我觉得一定要和不定积分联系起来，因为很多方法都和不定积分相似，比如说换元法，分部法等。

Integrate[Sqrt[r^2 - x^2]/(x + a), {x, -r, r},

Assumptions ->{0 <r <1, a >10 r}]

这是积分

Plot3D[\[Pi] (a - Sqrt[a^2 - r^2]), {r, -8, 8}, {a, -8, 8}]

Manipulate[Plot[\[Pi] (a - Sqrt[a^2 - r^2]), {r, -8, 8}], {a, -8, 8}]

这两个是画图

题主的做法主要存在两点问题：

1、quad函数用于计算数值积分，函数表达式中不能包含符号量；

2、被积函数的表达式应该写成关于被积变量的向量化的形式（也就是应该用点运算）。

参考代码：

R=1

syms L

rr = 0 : 0.1 : 1

for ii = 1 : length(rr)

    r = rr(ii)

    f = @(l)(acos((1+l*l-r*r)/(2*l))+r*r*acos((r*r+l*l-1)/(2*r*l))-0.5*sqrt(4*r*r-(1+r*r-l*l)^2))*2*l/(pi*r^4)

    fun = @(L) arrayfun(f,L)

    J(ii) = quadl(fun,0,r)

end

plot(rr, J)

或者也可以借用楼上 枫箫1 的部分代码，写成：

R=1

syms L

rr = 0 : 0.1 : 1

for ii = 1 : length(rr)

    r = rr(ii)

    SOA=R^2*acos((R^2+L^2-r^2)/(2*R*L))+r^2*acos((r^2+L^2-R^2)/(2*r*L))-...

        0.5*sqrt(4*R^2*r^2-(R^2+r^2-L^2)^2)

    PAB=SOA/(pi*r^2)

    p=2*L/r^2

    f=PAB*p

    fun = eval(['@(L)' vectorize(f)])

    fun = @(l) arrayfun(@(L)eval(f),l)

    J(ii) = quadl(fun,0,r)

end

plot(rr, J)

以上代码虽可以运行，但被积函数存在问题——SOA的第一项反余弦的值可能出现复数（因为在r稍小的情况下，acos的参数大于1），请题主再仔细检查一下。