Spectral Method伪代码实现Legendre多项式及其导数的生成

#include<iostream>
#include <cmath>
using namespace std;
double legendre(int N,int M,double x){
	//cout<<"leg"<<endl;
	double r;
	double s[10][3];
	int k=0,j=0;
	for(j=0;j<=M;j++){
		if(j==0){
			s[0][j]=1;
			s[1][j]=x;
			//cout<<s[1][j]<<endl;
			for(int k=1;k<=N-1;k++){
				s[k+1][j]=((2*k+1)*x*s[k][j]-k*s[k-1][j])/(k+1);
				//cout<<s[k+1][j]<<endl;
			}
		}else{
			s[0][j]=0;
		    if(j==1)
				s[1][j]=1;
			else
				s[1][j]=0;
			for(k=1;k<=N-1;k++){
			s[k+1][j]=(2*k+1)*s[k][j-1]+s[k-1][j];}
		}

	}
	r=s[N][M];
	//cout<<r<<"r"<<endl;
	return r;
}

double *lguass1(double p[10],int N,int M,double c){
	//cout<<"lguass1"<<endl;
	double z[10];
	double H=(double)1/49,a=-1.0,b=0.0;
	double x=0.0,xright=0.0;

	for(int k=1;k<=N-M;k++){
		b=a+H;
		//cout<<a<<b<<endl;
		while((legendre(N,M,a)*legendre(N,M,b))>0){
			a=b;
			b=a+H;
		}
		x=(a+b)/2;
		xright=b;
		while(fabs(x-xright)>=c){
			xright=x;
			x=x-legendre(N,M,x)/legendre(N,M+1,x);
			//cout<<2<<endl;
		}
		z[k-1]=x;
		a=x+H;
	    //cout<<x<<endl;

	}
       for(int i = 0;i<=N-M-1; i++)
       p[i] = z[i];
       return p;
}


int main(){
	int N;
	int M;
	double c;
    double *p = new double[10];	
	cout<<"N is the number of Legendre polynomial";
	cin>>N;
	cout<<"M is the order of derivative";
	cin>>M;
	cout<<"c is the accuracy tolerence";
	cin>>c;

    p =lguass1(p,N,M,c);
    for(int i=0;i<=N-M-1;i++){
		cout << p[i] << endl;
	}
    
    delete p;
	return 0;
} 

function z=guassmain
N = input('Enter N: ');
m = input('Enter m is the order of derivative: ');
c = input('Enter c is the accuracy tolerence : ');
%N=7;
%m=0;
%c=1e-8;
z=lgauss1(N,c,m);
end


function z=lgauss1(N,c,m)%the zeros of L_{N}^{(m)}(x)
%zeros of Legendre polynomials
%input N,c,m
%c is the accutacy tolerence
H=N^(-2);
a=-1;
for k=1:N-m
    %the following is to search the small interval containing a root
    b=a+H;
    while Legendre(N,m,a).*Legendre(N,m,b)>0
        a=b;b=a+H;
    end
    %the Newton's method in (a,b)
    x=(a+b)./2;
    xright=b;
   while abs(x-xright)>c||abs(x-xright)==c
        xright=x;
        x=x-Legendre(N,m,x)./Legendre(N,m+1,x);
   end
    z(k)=x;
    a=x+H;%move to anther interval containing a root
end
end


function r=Legendre(n,m,x)
%(n,m)：n次legendre函数的m次导数值
s=zeros(n+1,m+1);
for j=1:m+1
    if j==1
        s(1,j)=1;s(2,j)=x;%s(0,0)=1,s(1,0)=x
        for k=2:n
            s(k+1,j)=((2*(k-1)+1).*x.*s(k,j)-(k-1).*s(k-1,j))./((k-1)+1);
            %(n+1)s(k,0)=(2n+1)xs(n,0)-ns(n-1,0) (1.3.18)
        end
    else
        s(1,j)=0;%s(0,1)=s(0,2)=...=s(0,m)=0;
        if j==2
            s(2,j)=1;%s(1,1)=1
        else
            s(2,j)=0;%s(1,2)=s(1,3)=...=s(1,m)=0;
        end
        for k=2:n
            s(k+1,j)=(2.*(k-1)+1).*s(k,j-1)+s(k-1,j);
            %L_{n}(x)=1/(2n+1)(L'_{n+1}(x)-L'_{n-1}(x)) (1.3.22b)
        end
    end
end
r=s(n+1,m+1);
end