我正在研究一种可以使用基本线段计算曲线长度的算法。所以若我有曲线的x和y坐标的向量,我需要计算这个基本线段的数量。为此,我发明了自己的递归算法。有一个代码:
#include <iostream>
#include <vector>
#include <cmath>
using namespace std;
//***********************************************
inline double Pitagoras(double xp, double yp, double xf, double yf)
{
return sqrt((xp - xf)*(xp - xf) + (yp - yf)*(yp - yf));
}
//***************************************************
inline double calculateX(double xp, double yp, double a, double b, const double eps)
{
double delta;
double x1, x2;
delta = (-2.0 * xp + 2.0 *a*b - 2.0 *a*yp)*(-2.0 * xp + 2.0 *a*b - 2.0 *a*yp)
- 4.0* (1.0 + a*a)*(xp*xp + yp*yp + b*b - 2.0*b*yp - eps*eps);
x1 = (-(-2.0 * xp + 2.0 *a*b - 2.0 *a*yp) - sqrt(delta))/(2.0 * (1.0 + a*a));
x2 = (-(-2.0 * xp + 2.0 *a*b - 2.0 *a*yp) + sqrt(delta))/(2.0 * (1.0 + a*a));
if(x1 >= xp)
return x1;
else
return x2;
}
//***************************************************
inline double calculateY(double x, double a, double b)
{
return a*x + b;
}
//***********************************************
unsigned long algorithmKolmogorow(double xp, double yp, double xf,
double yf, const double eps, vector<double> &vectorX, vector<double> &vectorY);
//***********************************************
int main()
{
vector<double> vctrY; //vector of value of function
vector<double> vctrX; //vector of x
double xP,yP,xF,yF; //coordinates of two points on the curve
const double Eps = 0.0001; //length of elementary line
for(double x=1.0; x<=5 ;x +=0.001)
{
vctrX.push_back(x);
vctrY.push_back(x*x); //f(x) = x^2
}
xP = vctrX[0];
yP = vctrY[0];
xF = vctrX[1];
yF = vctrY[1]; //set beginning value
cout<<algorithmKolmogorow(xP, yP, xF, yF, Eps, vctrX, vctrY)*Eps;
return 0;
}
//***************************************************
unsigned long algorithmKolmogorow(double xp, double yp, double xf,
double yf, const double eps, vector<double> &vectorX, vector<double> &vectorY)
{
static unsigned long N; //licznik
static unsigned long i = 1;
double d;
double a,b;
d = Pitagoras(xp, yp, xf, yf);
if(d >= eps){
a = (yf - yp)/(xf - xp);
b = yp - a*xp;
xp = calculateX(xp, yp, a, b, eps);
yp = calculateY(xp, a, b);
N++;
}
else{
i++;
xf = vectorX[i];
yf = vectorY[i];
//cout<<i<<"t"<<vectorX[i]<<"t"<<vectorY[i]<<endl;
}
if(i < vectorX.size())
N = algorithmKolmogorow(xp, yp, xf, yf, eps, vectorX, vectorY);
return N;
}
正如你所看到的,主要是,我正在为抛物函数设置x,y坐标。当Eps很大时,例如Eps=0.001,一切都正常。如果我设置了一个较小的值,比如Eps=0.0001,那么我就会出现主题中的错误,程序就会停止运行。我完全不知道为什么。
我可以添加任何你需要的新信息(关于我的编译器,IDE,操作系统等(。
阵风
由于递归,您可能会耗尽堆栈空间。
一个快速的破解方法是在外壳中使用ulimit
来提高堆栈大小,例如
(ulimit -s unlimited; ./my_program)
然而,真正的解决方案是删除递归。algorithmKolmogorow
看起来只执行尾部递归,它总是可以转换为循环:
unsigned long algorithmKolmogorow(double xp, double yp, double xf,
double yf, const double eps, vector<double> &vectorX, vector<double> &vectorY)
{
static unsigned long N; //licznik
static unsigned long i = 1;
double d;
double a,b;
while(true) {
d = Pitagoras(xp, yp, xf, yf);
if(d >= eps){
a = (yf - yp)/(xf - xp);
b = yp - a*xp;
xp = calculateX(xp, yp, a, b, eps);
yp = calculateY(xp, a, b);
N++;
}
else{
i++;
if(i >= vectorX.size())
return N;
xf = vectorX[i];
yf = vectorY[i];
//cout<<i<<"t"<<vectorX[i]<<"t"<<vectorY[i]<<endl;
}
}
}
它还解决了访问一个元素通过其边界的向量的问题。
仍然存在一些问题:
- 应初始化
N
- 每次程序调用只能使用
algorithmKolmogorow
一次。从N
和i
中删除static
以修复此问题
i++;
xf = vectorX[i];
yf = vectorY[i];
//cout<<i<<"t"<<vectorX[i]<<"t"<<vectorY[i]<<endl;
}
if(i < vectorX.size())
这将超过向量的边界1。也许你需要if(i < vectorX.size()-1)
,或者(注意AndyG的评论(也许你需要将i++
放在i
的这两个用法之下。
Eps
的值太小。它会导致算法进展太慢,这意味着在计算完成之前,堆栈空间就用完了。我将Eps
增加到0.005,代码运行良好。但不确定它是否产生了正确的答案;你得检查一下。
此外,您需要更改这些行:
xf = vectorX[i];
yf = vectorY[i];
对此:
xf = vectorX[i-1];
yf = vectorY[i-1];
否则,您将访问上一次迭代中的最后一个向量元素。您可能还想对Eps
进行公差研究,看看它必须有多小才能得到准确的计算。