python解决缔合勒让德函数的递推(标准向前行递推算法)
在求解拉普拉斯面球谐函数时存在勒让德方程,而勒让德函数为勒让德方程的解。勒让德函数,通常以符号表 示,n为阶数,它满足于勒让德微分方程(1)。
完全规格化缔合勒让德函数是对缔合勒让德函数进行规格化后的结果,而则是通过勒让德函数计算得到的。θ为球坐标系下的一个坐标量,θ∈[0°,180°],也称为天顶距;n、m均为0或正整数,分别称为函数的阶、次,且n≥m。常用的递推算法主要有4种:标准向前列递推算法、标准向前行递推算法、跨阶次递推算法和Belikov递推算法。
在本次作业中,我使用的是“标准向前行递推算法”,其算法模型如下式:
式中,n≥2且n>m,其余参数如下:
采用python编写,输入所需的n和m,计算θ=36°时对应的,迭代至n=1000,代码如下:
1 def Lengendre(): 2 n = float(input('n=')) 3 m = float(input('m=')) 4 theta = float(36) 5 theta = theta*math.pi/180 6 p = np.zeros((1005, 1005)) 7 p[0][0] = float(1) 8 p[1][0] = math.sqrt(3)*math.cos(theta) 9 p[1][1] = math.sqrt(3)*math.sin(theta) 10 for i in range(2,1001): 11 for j in range(i,-1,-1): 12 if i == j: 13 c_mm = math.sqrt((2*j+1)/(2*j)) 14 p[i][i] = c_mm * math.sin(theta)*p[i-1][i-1] 15 else: 16 if j == 0: 17 j_j = float(2) 18 else: 19 j_j = float(1) 20 h_nm = math.sqrt(((i+j+2)*(i-j-1))/((i-j)*(i+j+1))) 21 g_nm = 2*(j+1)/(math.sqrt((i-j)*(i+j+1))) 22 if h_nm == 0: 23 p[i][j]=(1/math.sqrt(j_j)) * g_nm * math.cos(theta)/math.sin(theta)*p[i][i] 24 else: 25 p[i][j]=(1/math.sqrt(j_j)) * (g_nm * math.cos(theta)/(math.sin(theta))*p[i][j+1] - (h_nm * p[i][j+2])) 26 27 print(' L ',' m ',' P ') 28 print(int(n),int(m),p[int(n)][int(m)])