我正在尝试编写一个应用数学程序,用于计算复平面上的等高线积分。首先,我想编写梯形法的算法,但我对其外观有些困惑。毕竟 - 我们通常认为梯形法是针对二维图像的,而这里我们有f:C->C,所以我们在谈论四维。
最终,我希望使用此算法计算残留值,但当我尝试最简单的f(z)= 1 / z并将轮廓线设置为围绕原点的半径为1的圆时,我得不到接近1的值(这是我应该得到的)。以下是我的梯形法代码:
最终,我希望使用此算法计算残留值,但当我尝试最简单的f(z)= 1 / z并将轮廓线设置为围绕原点的半径为1的圆时,我得不到接近1的值(这是我应该得到的)。以下是我的梯形法代码:
complexCartesian trapezoid(complexCartesian *c1, complexCartesian *c2)
{
complexCartesian difference = *c1 - *c2;
complexCartesian ans(0.5 * (function(c1).real + function(c2).real) *
difference.mag(),
0.5 * (function(c1).imag + function(c2).imag) *
difference.mag());
return ans;
}
这里的 'function' 函数计算的是 f(z) = 1/z (我确定计算正确),而 complexCartesian 则是用于表示 a + bi 格式下的复数点的类:
class complexCartesian
{
public:
double real;
double imag;
complexCartesian operator+ (const complexCartesian& c) const;
complexCartesian operator- (const complexCartesian& c) const;
complexCartesian(double realLocal, double imagLocal);
double mag(); //magnitude
string toString();
complexPolar toPolar();
};
我很有信心这些函数都能够正常运行。(我知道有一个标准的复数类,但为了练习,我决定自己写。) 要进行积分,我使用以下方法:
double increment = .00001;
double radius = 1.0;
complexCartesian res(0,0); //result
complexCartesian previous(1, 0); //start the point 1 + 0i
for (double i = increment; i <= 2*PI; i+=increment)
{
counter++;
complex cur(radius * cos(i), radius * sin(i));
res = res + trapezoid(&cur, &previous);
previous = cur;
}
cout << "computed at " << counter << " values " << endl;
cout << "the integral evaluates to " << res.toString() << endl;
当我只沿着实轴进行积分,或者将我的函数替换为常数时,我得到正确的结果。否则,我往往会得到10^(-10)到10^(-15)数量级的数字。如果您有任何建议,我会非常感激。
谢谢。