勒让德多项式计算

简单写一下这次物理大地的程序吧，其实程序比较简单，主要还是要弄懂用哪个公式。
这是这次的作业要求：

计算勒让德公式可以用下面这个公式，当然是算不出来的，因为数字太大溢出了：
所以，我们用下面这个公式进行迭代：

具体代码如下：
下面是用python写的

# -*- coding: utf-8 -*-
import math


def claulate_legendre(n, sita_rad):
    P = list()
    for row in range(n+1):
        P.append([])
        for clumn in range(row+1):
            P[row].append(0)
    P[0][0] = 1
    P[1][0] = math.sqrt(3) * math.cos(sita_rad)
    P[1][1] = math.sqrt(3) * math.sin(sita_rad)
    for row in range(2, n+1):
        coefficient_1 = math.sqrt((2 * row + 1)/(2 * row))
        coefficient_2 = math.sqrt(2 * row + 1)
        P[row][row] = coefficient_1 * math.sin(sita_rad) * P[row-1][row-1]
        P[row][row-1] = coefficient_2 * math.cos(sita_rad) * P[row-1][row-1]
    for row in range(2, n+1):
        for clumn in range(0, row-1):
            coefficient_3 = math.sqrt((2 * row + 1) / (row - clumn) / (row + clumn))
            coefficient_4 = math.sqrt(2 * row - 1)
            coefficient_5 = math.sqrt((row - clumn - 1) * (row + clumn -1) / (2 * row - 3))
            P[row][clumn] = coefficient_3 * (coefficient_4 * math.cos(sita_rad) * P[row - 1][clumn]
                                             - coefficient_5 * P[row - 2][clumn])
    return P


def write_file(out_filename, P_data, n):
    lines = list()
    header_line = "degree order value\n"
    lines.append(header_line)
    for m in range(0, n+1):
        line = "{:<7d}{:<6d}{:.16f}\n".format(n, m, P_data[n][m])
        lines.append(line)
    with open(out_filename, 'w') as fopen:
        fopen.writelines(lines)


def main():
    id_number = 293
    n = id_number
    sita = [75, 25]
    out_filename = ["C:\\Users\\王阳\\Desktop\\2017301610293_75.txt", "C:\\Users\\王阳\\Desktop\\2017301610293_25.txt"]
    for i in range(2):
        sita_rad = math.radians(sita[i])
        P_data = claulate_legendre(n, sita_rad)
        write_file(out_filename[i], P_data, n)


if __name__ == "__main__":
    main()

下面是用C++写的，结果是一样的

#include<iostream>
#include<fstream>
#include<iomanip>
using namespace std;
#define PI 4*atan(1.0)

void test(int n,double sita,string filename)
{
	double sita_rad=sita*PI/180;
	double **P=new double* [n+1];
	for(int i=0;i<n+1;i++)
	{
		P[i]=new double[n+1];
	}
	P[0][0]=1;
	P[1][0]=sqrt(3.0)*cos(sita_rad);
	P[1][1]=sqrt(3.0)*sin(sita_rad);
	for(int i=2;i<=n;i++)
	{
		P[i][i]=sqrt((2*i+1.0)/(2*i))*sin(sita_rad)*P[i-1][i-1];
		P[i][i-1]=sqrt(2*i+1.0)*cos(sita_rad)*P[i-1][i-1];
	}
	for(int i=2;i<=n;i++)
	{
		for(int j=0;j<=i-2;j++)
		{
			double f3=sqrt((2*i+1.0)/(i-j)/(i+j));
			double f4=sqrt(2*i-1.0);
			double f5=sqrt((i-j-1)*(i+j-1)/(2*i-3.0));
			P[i][j]=f3*(f4*cos(sita_rad)*P[i-1][j]-f5*P[i-2][j]);
		}
	}

	ofstream ofile;
	ofile.open(filename);
	ofile<<"degree order value"<<endl;
	for(int i=0;i<n+1;i++)
	{
		ofile<<setprecision(16)<<setiosflags(ios::left);
		ofile<<setw(7)<<n<<setw(6)<<i<<setw(16)<<P[n][i]<<endl;
	}
	ofile.close();
}

int main(void)
{
	string filename="C:\\Users\\王阳\\Desktop\\2017301610293_75.txt";
	int id_number=293;
	int n=id_number;
	double sita = 75;
	test(n,sita,filename);
	system("pause");
	return 0;
}