各位虫友:
我遇到了一个比较复杂的2重积分,我的积分函数表达式如下:
如图片所示
在对上面的kx及ky求二重积分时无法用mathematica中的语句直接计算,因为积分函数相当复杂,用直接的mathematica语句无法积分,在网上查阅了一个用c++编写出来的版本,现在我想把他转换成mathematica的语句。希望对此熟悉的朋友能帮忙转换一下。
c++语句如下:
/* ============================================================= *
* trapzd2d. Computes the nth stage of refinement of an extended *
* 2d trapezian rule. func is a pointer to a function with two *
* real variables, it is to be integrated *
* between ax<=x<=ay and ay<=y<=by. Works analogous to the one- *
* dimensional trapzd (Numerical Recipes, p.137). *
* ============================================================= */
typedef double real; /* can be put e.g. to float if single precission is sufficient */
#define FUNC(x,y) ((*func)(x,y))
real trapzd2d(real (*func)(real, real), real ax, real bx, real ay, real by, int n)
{
real x, y, tnm, sum, delx, dely;
static real s;
int it, j, k;
if (n==1) {
/* 1. stage: Function is evaluated only at the corners. */
return (s=0.25*(bx-ax)*(by-ay)*(FUNC(ax,ay)+FUNC(bx,ay)+FUNC(ax,by)+FUNC(bx,by)));
} else {
for (it=1, j=1; j<n-1; j++) it <<= 1;
tnm = it;
/* delx, dely are the spacings of the points to be added */
delx = (bx-ax)/tnm;
dely = (by-ay)/tnm;
x = ax + 0.5*delx;
/* Evaluates function at the borders of the integration area */
for (sum=0.0, j=1; j<=it; j++, x+=delx) {
sum += 0.5*(FUNC(x,ay) + FUNC(x,by));
}
y = ay + 0.5*dely;
for (j=1; j<=it; j++, y+=dely) {
sum += 0.5*(FUNC(ax,y) + FUNC(bx,y));
}
/* Interiour of the area */
x = ax + 0.5*delx;
for (j=1; j<=it; j++, x+= delx) {
y = ay + 0.5*dely;
for (k=1; k<=it; k++, y+=dely) {
sum += FUNC(x,y);
}
}
x = ax + delx;
for (j=1; j<it; j++, x+=delx) {
y = ay + 0.5*dely;
for (k=1; k<=it; k++, y+=dely) {
sum += FUNC(x,y);
}
}
x = ax + 0.5*delx;
for (j=1; j<=it; j++, x+=delx) {
y = ay + dely;
for (k=1; k<it; k++, y+=dely) {
sum += FUNC(x,y);
}
}
/* Replace s by its refined value */
s = 0.25*(s + (bx-ax)*(by-ay)*sum/(tnm*tnm));
return s;
}
}
#undef FUNC
或者各位虫友能有更好的解决这个积分的办法,也可以麻烦告知,多谢
![mathematica中如何实现c++程序转换为mathematica语句]()
积分函数.png |