C++：最速降线问题（用两个已知点求待定系数，并输出）

最新推荐文章于 2023-03-18 15:23:25 发布

水元素法师

最新推荐文章于 2023-03-18 15:23:25 发布

阅读量489

点赞数 2

分类专栏：兴趣使然的水文文章标签：二分法 windows 抽象代数

本文链接：https://blog.csdn.net/weixin_44960715/article/details/106016929

版权

兴趣使然的水文专栏收录该内容

8 篇文章 0 订阅

订阅专栏

其实没啥内容好说明的，大可以只看源码，下载编译运行，可能更直观。
windows 或 Linux 下编译皆可，无数学库安装要求。

源码链接，0积分下载

#include <iostream>
#include <cmath>
#define g 9.8
#define step 0.1
#define accuracyForDichotomy 1e-6
enum
{
	section = 5
}; //In one direction
struct point
{
	double x;
	double y;
};
double xx0, yy0, a;						   //Undetermined Coefficient
double nullPointAbscissa[2 * section + 2]; //10 section: -5,-4,...,4,5
double extremum[2 * section];

//(cos(x)-1)/(x-sin(x))=ratio
//2cos(x)+x*sin(x)=2
void getNullPointAbscissa(void); //prepare for getUndeterminedCoefficient()

//x=xx0+a(t-sin(t)); y=yy0+a(cos(t)-1)
double f(double);

//ff0(x)=(cos(x)-1)/(x-sin(x))
double ff0(double x) { return (cos(x) - 1) / (x - sin(x)); };

//ff1(x)=x-sin(x)
double ff1(double x) { return x - sin(x); };

//ff2(x)=2cos(x)+xsin(x)-2
double ff2(double x) { return 2 * cos(x) + x * sin(x) - 2; };

//dichotomy
double dichotomy(double x1, double x2, double value, double (*pf)(double));
void getUndeterminedCoefficient(point, point);

using namespace std;
int main(void)
{
	FILE *outPut = fopen("Curve.txt", "w");
	point p1, p2;
	double xxx;

	cout << "p1.x=";
	cin >> p1.x;
	cout << "p1.y=";
	cin >> p1.y;
	cout << "p2.x=";
	cin >> p2.x;
	cout << "p2.y=";
	cin >> p2.y;
	cout << endl;

	getNullPointAbscissa();										 //Prepare for getUndeterminedCoefficient. Set some arrays.
	getUndeterminedCoefficient(p1, p2);							 //set xx0,yy0,a
	cout << "x0=" << xx0 << " y0=" << yy0 << " a=" << a << endl; //show xx0,yy0,a

	xxx = p1.x;
	while (xxx < p2.x)
	{
		fprintf(outPut, "%f %f\n", xxx, f(xxx));
		xxx += step;
	}

	fclose(outPut);
	//for(int i=0;i<2*section+2;i++ ) cout<<nullPointAbscissa[i]<<" ";
	//cout << endl;
	//for(int i=0;i<2*section;i++ ) cout<<extremum[i]<<" ";
	return 0;
}

void getNullPointAbscissa(void)
{
	const double interval = 1e-3;
	double xx;
	int count = 0;
	bool last, next; //True->up False->down
	//Forward traversal
	xx = interval;
	last = ff2(xx) > 0;
	while (count < section)
	{
		xx += interval;
		next = ff2(xx) > 0;
		if (last != next) //Different side.
		{
			count++;
			nullPointAbscissa[count + section + 1] = xx - interval / 2;
			extremum[count + section - 1] = ff0(nullPointAbscissa[count + section + 1]);
		}
		last = next;
	}
	//Function is symmetrical
	nullPointAbscissa[section] = -0.1;
	nullPointAbscissa[section + 1] = 0.1;
	for (int i = 0; i < section; i++)
	{
		nullPointAbscissa[i] = -nullPointAbscissa[2 * section + 1 - i];
		extremum[i] = ff0(nullPointAbscissa[i]);
	}
}

void getUndeterminedCoefficient(point p1, point p2)
{
	double ratio, tempV, theta;
	int whichArea;

	xx0 = p1.x;
	yy0 = p1.y;

	//Get coefficient a
	ratio = (p2.y - yy0) / (p2.x - xx0);

	//According to the value of ratio, give the user optional interval mark number.
	if (ratio > 1e-15)
	{
		cout << "optional mark number of the interval: -1, ";
		for (int i = section - 2; i >= 0; i--)
		{
			if ((ratio - extremum[i]) * (ratio - extremum[i + 1]) < 0) //Between
				cout << ", " << i - section;
		}
		cout << endl
			 << "And you have just chosen: ";
		cin >> whichArea;
	}
	else if (ratio < -1e-15)
	{
		cout << "optional mark number of the interval: 1";
		for (int i = section; i < 2 * section - 1; i++)
		{
			if ((ratio - extremum[i]) * (ratio - extremum[i + 1]) < 0) //Between
				cout << ", " << i - section + 2;
		}
		cout << endl
			 << "And you have just chosen: ";
		cin >> whichArea;
	}
	else //ratio -> 0
	{
		const double Pi = 3.141592653589793;
		cout << "(y2-y1)/(x2-x1) is too small, just enter 1,2,... for example:";
		cin >> whichArea;
		theta = 2 * Pi * whichArea;
		a = ((p2.x - xx0) / (theta - sin(theta)) + (p2.y - yy0) / (cos(theta) - 1)) / 2;
		return;
	}
	//Using dichotomy to get theta
	if (whichArea > 0)
	{
		theta = dichotomy(nullPointAbscissa[whichArea + section], nullPointAbscissa[whichArea + section + 1], ratio, ff0);
	}
	else
	{
		theta = dichotomy(nullPointAbscissa[whichArea + section], nullPointAbscissa[whichArea + section + 1], ratio, ff0);
	}
	//Set a
	a = ((p2.x - xx0) / (theta - sin(theta)) + (p2.y - yy0) / (cos(theta) - 1)) / 2;
}

double f(double x)
{
	double theta, value;

	value = (x - xx0) / a;
	theta = dichotomy(value - 1.1, value + 1.1, value, ff1);

	return yy0 + a * (cos(theta) - 1);
}

double dichotomy(double x1, double x2, double value, double (*pf)(double))
{
	if ((value - (*pf)(x1)) * (value - (*pf)(x2)) > 0)
	{
		cout << "ERROR: In dichotomy, value conflicting.\n";
		return -1;
	}
	double xm;
	bool leftB, rightB;

	do
	{
		xm = (x1 + x2) / 2;

		leftB = (value - (*pf)(x1)) * (value - (*pf)(xm)) < 0;	//value is in the left interval
		rightB = (value - (*pf)(xm)) * (value - (*pf)(x2)) < 0; //value is in the right interval

		if (leftB && rightB)
		{
			cout << "ERROR: In dichotomy, Function is non-monotonic.\nYou're not that straight! LOL!\n";
			return -1;
		}
		else if (leftB || rightB) //If running to here, it means 0 1 or 1 0 or 0 0.
		{
			//cout << "Dichotomy: Temporarily Normal\n";
			if (leftB)
			{

				x2 = xm;
			}
			else
			{
				x1 = xm;
			}
		}
		else //If running to here, it can only be 0 0.
		{
			string str;
			cout << "ERROR: In dichotomy, a rarely weird error pop up.\nDid you let a=b? [Y/N] ";
			cin >> str;
			if (str == "Y" || str == "y")
				return 0; //This 0 is for f(p1.x)
			else
				return -1;
		}

	} while (fabs((*pf)(xm)-value) > accuracyForDichotomy);

	return xm;
}

说明文档

1 算法流程

1.1 主函数

int main(void)
{
	FILE *outPut = fopen("Curve.txt", "w");
	point p1, p2;
	double xxx;

	cout << "p1.x=";
	cin >> p1.x;
	cout << "p1.y=";
	cin >> p1.y;
	cout << "p2.x=";
	cin >> p2.x;
	cout << "p2.y=";
	cin >> p2.y;
	cout << endl;

	getNullPointAbscissa();										 
    //Prepare for getUndeterminedCoefficient. Set some arrays.
	getUndeterminedCoefficient(p1, p2);							 
    //set xx0,yy0,a
	cout << "x0=" << xx0 << " y0=" << yy0 << " a=" << a << endl; 
    //show xx0,yy0,a

	xxx = p1.x;
	while (xxx < p2.x)
	{
		fprintf(outPut, "%f %f\n", xxx, f(xxx));
		xxx += step;
	}

	fclose(outPut);
	//for(int i=0;i<2*section+2;i++ ) cout<<nullPointAbscissa[i]<<" ";
	//cout << endl;
	//for(int i=0;i<2*section;i++ ) cout<<extremum[i]<<" ";
	return 0;
}

$\begin{cases} x=x_0+a(\theta-sin\theta) \\ y=y_0+a(cos\theta-1) \end{cases}$

平摆线参数方程如上，其中 $x_0,y_0$ 即为第一个点的横纵坐标。用另一个坐标求出a的值。

1.2 副函数

//(cos(x)-1)/(x-sin(x))=ratio
//2cos(x)+x*sin(x)=2
void getNullPointAbscissa(void); //prepare for getUndeterminedCoefficient()

//x=xx0+a(t-sin(t)); y=yy0+a(cos(t)-1)
double f(double);

//ff0(x)=(cos(x)-1)/(x-sin(x))
double ff0(double x) { return (cos(x) - 1) / (x - sin(x)); };

//ff1(x)=x-sin(x)
double ff1(double x) { return x - sin(x); };

//ff2(x)=2cos(x)+xsin(x)-2
double ff2(double x) { return 2 * cos(x) + x * sin(x) - 2; };

//dichotomy
double dichotomy(double x1, double x2, double value, double (*pf)(double));
void getUndeterminedCoefficient(point, point);

1.2.1 参数的含义

#define g 9.8
#define step 0.1
#define accuracyForDichotomy 1e-6
enum
{
	section = 5
}; //In one direction
struct point
{
	double x;
	double y;
};
double xx0, yy0, a;						   //Undetermined Coefficient
double nullPointAbscissa[2 * section + 2]; //10 section: -5,-4,...,4,5
double extremum[2 * section];

g 为重力加速度。
step 为输出f(x)的x的步长。
accuracyForDichotomy 为二分法的精度要求。
section 为所取的正向或负向的单调域（函数 ff0(x) ）个数。
结构 point 更方便书写。
xx0, yy0, a 为参数方程里的待定系数。
数组 nullPointAbscissa[2 * section + 2] 存储单调域端点的横坐标，零附近取-0.1和0.1两个值，故要求ratio < 30，这个条件一般都是符合的。
数组 extremum[2 * section] 存储非零端点的函数值。

1.2.2 函数ff0(x), ff1(x)和ff2(x)的图像

//ff0(x)=(cos(x)-1)/(x-sin(x))
double ff0(double x) { return (cos(x) - 1) / (x - sin(x)); };

//ff1(x)=x-sin(x)
double ff1(double x) { return x - sin(x); };

//ff2(x)=2cos(x)+xsin(x)-2
double ff2(double x) { return 2 * cos(x) + x * sin(x) - 2; };

颜色分别为：深蓝、棕褐、荧光绿。

p.s. 青蓝色是 $\frac{d}{dx}ff0(x)$
嘿嘿，就是这么简陋
难道还有比这更简陋的作图方式吗

ff2是ff0的导函数的分子，令其为零，得极值点。

1.2.3 核心副函数

void getNullPointAbscissa(void); //prepare for getUndeterminedCoefficient()
void getUndeterminedCoefficient(point, point);
double dichotomy(double x1, double x2, double value, double (*pf)(double));
double f(double);                //x=xx0+a(t-sin(t)); y=yy0+a(cos(t)-1)

1.2.3.1 获取导函数零点坐标getNullPointAbscissa

找的是 ff2(x) 绿色函数线的零点。考虑到单调区间不便于划分，直接遍历找变号零点。

void getNullPointAbscissa(void)
{
	const double interval = 1e-3;
	double xx;
	int count = 0;
	bool last, next; //True->up False->down
	//Forward traversal
	xx = interval;
	last = ff2(xx) > 0;
	while (count < section)
	{
		xx += interval;
		next = ff2(xx) > 0;
		if (last != next) //Different side.
		{
			count++;
			nullPointAbscissa[count + section + 1] = xx - interval / 2;
			extremum[count + section - 1] = ff0(nullPointAbscissa[count + section + 1]);
		}
		last = next;
	}
	//Function is symmetrical
	nullPointAbscissa[section] = -0.1;
	nullPointAbscissa[section + 1] = 0.1;
	for (int i = 0; i < section; i++)
	{
		nullPointAbscissa[i] = -nullPointAbscissa[2 * section + 1 - i];
		extremum[i] = ff0(nullPointAbscissa[i]);
	}
}

1.2.3.2 getUndeterminedCoefficient：目的在于计算参数a

为了计算a，参照 (1) 式，可以定义ratio如下
$ratio=\frac{y_2-y_1}{x_2-x_1}=\frac{cos\theta-1}{\theta-sin\theta}$
注意到对于两个定点 $p_1 (x_1,y_1),p_2(x_2,y_2)，有$ $x_0=x_1, y_0=y_1$ 这一条件。

ratio 是一定值，主要是想求出 ratio 对应的 $\theta$ ，从而由
$a=\frac{x_2-x_1}{\theta-sin\theta} =\frac{y_2-y_1}{cos\theta-1} =\left(\frac{x_2-x_1}{\theta-sin\theta}+\frac{y_2-y_1}{cos\theta-1} \right)/2$
给出a的值。

关键在于求 $\theta$ ，本人采用的是二分法，在这里就会用到那两个数组。

void getUndeterminedCoefficient(point p1, point p2)
{
	double ratio, tempV, theta;
	int whichArea;

	xx0 = p1.x;
	yy0 = p1.y;

	//Get coefficient a
	ratio = (p2.y - yy0) / (p2.x - xx0);

	//According to the value of ratio, give the user optional interval mark number.
	if (ratio > 1e-15)
	{
		cout << "optional mark number of the interval: -1, ";
		for (int i = section - 2; i >= 0; i--)
		{
			if ((ratio - extremum[i]) * (ratio - extremum[i + 1]) < 0) //Between
				cout << ", " << i - section;
		}
		cout << endl
			 << "And you have just chosen: ";
		cin >> whichArea;
	}
	else if (ratio < -1e-15)
	{
		cout << "optional mark number of the interval: 1";
		for (int i = section; i < 2 * section - 1; i++)
		{
			if ((ratio - extremum[i]) * (ratio - extremum[i + 1]) < 0) //Between
				cout << ", " << i - section + 2;
		}
		cout << endl
			 << "And you have just chosen: ";
		cin >> whichArea;
	}
	else //ratio -> 0
	{
		const double Pi = 3.141592653589793;
		cout << "(y2-y1)/(x2-x1) is too small, just enter 1,2,... for example:";
		cin >> whichArea;
		theta = 2 * Pi * whichArea;
		a = ((p2.x - xx0) / (theta - sin(theta)) + (p2.y - yy0) / (cos(theta) - 1)) / 2;
		return;
	}
	
    //Using dichotomy to get theta
	if (whichArea > 0)
	{
		theta = dichotomy(nullPointAbscissa[whichArea + section], nullPointAbscissa[whichArea + section + 1], ratio, ff0);
	}
	else
	{
		theta = dichotomy(nullPointAbscissa[whichArea + section], nullPointAbscissa[whichArea + section + 1], ratio, ff0);
	}
	
    //Set a
	a = ((p2.x - xx0) / (theta - sin(theta)) + (p2.y - yy0) / (cos(theta) - 1)) / 2;
}

这里的判断都是根据函数特性的优化结果。

然后供选择的档位1, 2, …, section 是最终曲线的周期数，一般而言，whichArea=1 。

1.2.3.3 二分法 dichotomy

在这里的亮点是：

有报错补正功能。
运用函数指针，灵活度高。

double dichotomy(double x1, double x2, double value, double (*pf)(double))
{
	if ((value - (*pf)(x1)) * (value - (*pf)(x2)) > 0)
	{
		cout << "ERROR: In dichotomy, value conflicting.\n";
		return -1;
	}
	double xm;
	bool leftB, rightB;

	do
	{
		xm = (x1 + x2) / 2;

		leftB = (value - (*pf)(x1)) * (value - (*pf)(xm)) < 0;	
        //value is in the left interval
		rightB = (value - (*pf)(xm)) * (value - (*pf)(x2)) < 0; 
        //value is in the right interval

		if (leftB && rightB)
		{
			cout << "ERROR: In dichotomy, function is non-monotonic.\nYou're not that straight! LOL!\n";
			return -1;
		}
		else if (leftB || rightB) //If running to here, it means 0 1 or 1 0 or 0 0.
		{
			//cout << "Dichotomy: Temporarily Normal\n";
			if (leftB)
			{

				x2 = xm;
			}
			else
			{
				x1 = xm;
			}
		}
		else //If running to here, it can only be 0 0.
		{
			string str;
			cout << "ERROR: In dichotomy, a rarely weird error pop up.\nDid you let a=b? [Y/N] ";
			cin >> str;
			if (str == "Y" || str == "y")
				return 0; //This 0 is for f(p1.x)
			else
				return -1;
		}
	} while (fabs((*pf)(xm)-value) > accuracyForDichotomy);

	return xm;
}

1.2.3.4 参数方程转换：y=f(x)

double f(double x)
{
	double theta, value;

	value = (x - xx0) / a;
	theta = dichotomy(value - 1.1, value + 1.1, value, ff1);

	return yy0 + a * (cos(theta) - 1);
}

这里 value - 1.1, value + 1.1 是依据 ff1(x) = x-sin(x) 所作的优化。

思路还是二分求 $\theta$ ，回代 $\theta$ 得 y = f(x) 的值。

2 输出结果

对于 $p_1 (1,10), p_2(10,7)$ 两点，终端输入坐标：

win

输出 Curve.txt：

1.000000 10.000000
1.100000 9.573037
1.200000 9.332098
1.300000 9.135757
1.400000 8.965014
1.500000 8.811888
1.600000 8.672064
1.700000 8.542850
1.800000 8.422432
1.900000 8.309511
2.000000 8.203112
2.100000 8.102485
2.200000 8.007031
2.300000 7.916261
2.400000 7.829780
2.500000 7.747250
2.600000 7.668388
2.700000 7.592947
2.800000 7.520719
2.900000 7.451514
3.000000 7.385174
3.100000 7.321555
3.200000 7.260528
3.300000 7.201979
3.400000 7.145808
3.500000 7.091917
3.600000 7.040225
3.700000 6.990657
3.800000 6.943142
3.900000 6.897615
4.000000 6.854022
4.100000 6.812307
4.200000 6.772423
4.300000 6.734322
4.400000 6.697966
4.500000 6.663315
4.600000 6.630334
4.700000 6.598989
4.800000 6.569254
4.900000 6.541096
5.000000 6.514493
5.100000 6.489421
5.200000 6.465857
5.300000 6.443782
5.400000 6.423177
5.500000 6.404026
5.600000 6.386312
5.700000 6.370022
5.800000 6.355143
5.900000 6.341662
6.000000 6.329572
6.100000 6.318862
6.200000 6.309523
6.300000 6.301549
6.400000 6.294935
6.500000 6.289675
6.600000 6.285766
6.700000 6.283204
6.800000 6.281988
6.900000 6.282117
7.000000 6.283591
7.100000 6.286411
7.200000 6.290579
7.300000 6.296098
7.400000 6.302973
7.500000 6.311207
7.600000 6.320809
7.700000 6.331783
7.800000 6.344139
7.900000 6.357887
8.000000 6.373036
8.100000 6.389598
8.200000 6.407586
8.300000 6.427015
8.400000 6.447901
8.500000 6.470260
8.600000 6.494111
8.700000 6.519475
8.800000 6.546374
8.900000 6.574832
9.000000 6.604874
9.100000 6.636530
9.200000 6.669828
9.300000 6.704804
9.400000 6.741492
9.500000 6.779930
9.600000 6.820163
9.700000 6.862235
9.800000 6.906194
9.900000 6.952097
10.000000 7.000000

根据 Curve.txt 作图

粗制滥造

以上即为全部说明内容

水元素法师

关注

2
点赞
踩
5

收藏

觉得还不错? 一键收藏
0
评论
C++：最速降线问题（用两个已知点求待定系数，并输出）

其实没啥内容好说明的，大可以只看源码，下载编译运行，可能更直观。windows 或 Linux 下编译皆可，无数学库安装要求。源码链接，0积分下载#include <iostream>#include <cmath>#define g 9.8#define step 0.1#define accuracyForDichotomy 1e-6enum{ section = 5}; //In one directionstruct point{ double x
复制链接

扫一扫