缔合勒让德函数递推

最新推荐文章于 2024-09-13 16:27:15 发布

crazy小疯子

最新推荐文章于 2024-09-13 16:27:15 发布

阅读量912

点赞数 8

分类专栏： C++ 文章标签：算法 c++

本文链接：https://blog.csdn.net/qq_62772099/article/details/133249590

版权

C++ 专栏收录该内容

8 篇文章 4 订阅

订阅专栏

一：公式

标准向前列递推法

跨阶次递推法函数

代码

标准向前列递推法

void fun1(double o,int n)//标准向前列递推法，o为θ，n为要求的总阶数（设定可求的最大阶数为200）
//本次实习o=30度，n=181(0-180)
{
	double P[200][200];
	//种子点
	P[0][0] = 1;
	P[1][0] = sqrt(3)*cos(o);
	P[1][1] = sqrt(3)*sin(o);
	for (int l = 2; l < n; l++)
	{
		double f1 = sqrt(1 + 1.0 / 2 / l);
		double f2 = sqrt(2 * l + 1);//f1,f2与m无关，只与每一行的后两项递推有关
		for (int m = 0; m < l-1; m++)//每一行除了后两个的递推循坏，因为需要前面两行的值，这两列前面不够两个
		{
			double f3 = sqrt((2.0 * l + 1.0) / (l - m) / (l + m));
			double f4 = sqrt(2.0 * l - 1.0);
			double f5 = sqrt((l - m - 1) * (l + m - 1) *1.0/ (2.0* l - 3.0));//f3,f4,f5
			P[l][m] = f3 * (f4 * cos(o) * P[l - 1][m] - f5 * P[l - 2][m]);
		}

		P[l][l - 1] = f2 * cos(o) * P[l - 1][l - 1];//每一行倒数第二个的递推，用到这一列上一个的值
		P[l][l] = f1 * sin(o) * P[l - 1][l - 1];//每一行倒数第一个的递推，因为他是这一列的起点，递推用到前一行前一列的值
	}
	ofstream file;//文本输出流
	file.open("标准向前列递推法结果.txt");
	file << "l\t" << "m\t" << "Plm\t\n";
	for (int l = 0; l < n; l++)
	{
		for (int m = 0; m < l+1; m++)
		{
			//cout << l << "\t" << m << "\t" << setprecision(8) << P[l][m] << endl;
			file << l << "\t" << m << "\t" << setprecision(8) << P[l][m] << endl;//设置输出有效数位为5位
		}
	}

跨阶次递推法函数

void fun2(double o, int n)//跨阶次递推法，o为θ，n为要求的总阶数（设定可求的最大阶数为200）
//本次实习o=30度，n=181(0-180)
{
	double P[200][200];
	//种子点
	P[0][0] = 1;
	P[1][0] = sqrt(3) * cos(o);
	P[1][1] = sqrt(3) * sin(o);
	for (int l = 2; l < n; l++)
	{
		if (l == 2)//对于第三行，其第一列用前两行的值，第二列用上一行最后一个值(标准向前列递推法)
		{
			int m = 0; //其第一列用前两行的值
			double f3 = sqrt((2.0 * l + 1.0) / (l - m) / (l + m));
			double f4 = sqrt(2.0 * l - 1.0);
			double f5 = sqrt((l - m - 1.0) * (l + m - 1.0) * 1.0 / (2.0 * l - 3.0));
			P[l][m] = f3 * (f4 * cos(o) * P[l - 1][m] - f5 * P[l - 2][m]);
			m = 1;//第二列用上一行最后一个值,将 P[l - 2][m]设为0
			f3 = sqrt((2.0 * l + 1.0) / (l - m) / (l + m));
			f4 = sqrt(2.0 * l - 1.0);
			P[l][m] = f3 * f4 * cos(o) * P[l - 1][m];
			m = 2;//第三列用跨阶次递推法：系数要乘根号二，缺失的P[l - 2][m]设为0
			double b = sqrt(2.0 * (2.0 * l + 1.0) * (l + m - 2.0) * (l + m - 3.0) / (2.0 * l - 3.0) / (l + m) / (l + m - 1.0));
			double c = sqrt(2.0 * (l - m + 1.0) * (l - m + 2.0) / (l + m) / (l + m - 1.0));
			P[l][m] = b * P[l - 2][m - 2] - c * P[l][m - 2];
		}
		else//对于第四行起：前两列的值都是用前行的值(标准向前列递推法)
		{
			for (int m = 0; m < 2; m++)
		    {
				double f3 = sqrt((2.0 * l + 1.0) / (l - m) / (l + m));
		        double f4 = sqrt(2.0 * l - 1.0);
		        double f5 = sqrt((l - m - 1.0) * (l + m - 1.0) / (2.0 * l - 3.0));
			    P[l][m] = f3 * (f4 * cos(o) * P[l - 1][m] - f5 * P[l - 2][m]);
		    }
		    int m = 2;//对于m-2对于0的即第三列，使用跨阶次递推法，但是其bc系数要乘根号二。(跨阶次递推法)
			double a = sqrt((2.0 * l + 1.0) * (l - m) * (l - m - 1.0) / (2.0 * l - 3.0) / (l + m) / (l + m - 1.0));
			double b = sqrt(2.0 * (2.0 * l + 1.0) * (l + m - 2.0) * (l + m - 3.0) / (2.0 * l - 3.0) / (l + m) / (l + m - 1.0));
		    double c= sqrt(2.0 * (l - m + 1.0) * (l - m + 2.0)  / (l + m) / (l + m - 1.0));
	        P[l][m] = a * P[l - 2][m] + b * P[l - 2][m - 2] - c * P[l][m - 2];
	     	for (int m = 3; m < l - 1;m++)//每一行第四列起到倒数第二项，使用跨阶次递推法，其系数bc不用乘根号二(跨阶次递推法)
	      	{
			    double a = sqrt((2.0 * l + 1.0) * (l - m) * (l - m - 1.0) / (2.0 * l - 3.0) / (l + m) / (l + m - 1.0));
		     	double b = sqrt((2.0 * l + 1.0) * (l + m - 2.0) * (l + m - 3.0) / (2.0 * l - 3.0) / (l + m) / (l + m - 1.0));
		    	double c = sqrt((l - m + 1.0) * (l - m + 2.0) / (l + m) / (l + m - 1.0));
		    	P[l][m] = a * P[l - 2][m] + b * P[l - 2][m - 2] - c * P[l][m - 2];
		    }
			if (l == 3)//第四行，倒数第二个恰好为第三个数(m=2)前面已经求解，所以只剩下最后一个值
			{
				m = 3;
				double b = sqrt((2.0 * l + 1.0) * (l + m - 2.0) * (l + m - 3.0) / (2.0 * l - 3.0) / (l + m) / (l + m - 1.0));
				double c = sqrt((l - m + 1.0) * (l - m + 2.0) / (l + m) / (l + m - 1.0));
				P[l][m] = b * P[l - 2][m - 2] - c * P[l][m - 2];
			}
			else
			{
				for (int m = l - 1; m < l + 1; m++)//每一行最后两个值由于缺失其对应列前两行的值，所以将其值P[l - 2][m]设为0。
				{
					double b = sqrt((2.0 * l + 1.0) * (l + m - 2.0) * (l + m - 3.0) / (2.0 * l - 3.0) / (l + m) / (l + m - 1.0));
					double c = sqrt((l - m + 1.0) * (l - m + 2.0) / (l + m) / (l + m - 1.0));
					P[l][m] = b * P[l - 2][m - 2] - c * P[l][m - 2];
				}
			}
		}
	}
	ofstream file;//文本输出流
	file.open("跨阶次递推法结果.txt");
	file << "l\t" << "m\t" << "Plm\t\n";
	for (int l = 0; l < n; l++)
	{
		for (int m = 0; m < l + 1; m++)
		{
			//cout << l << "\t" << m << "\t" << setprecision(8) << P[l][m] << endl;
			file << l << "\t" << m << "\t" << setprecision(8) << P[l][m] << endl;//设置输出有效数位为5位
		}
	}