参考下面的程序吧
- % Determination of the coefficients in DQM
- % 19/03/2006 By Dr. Yang J., Modified by Mr. Xiang H.J.
- % Department of BC, City University of Hong Kong, Hong Kong
- % [C,X]=DQC(N,M,KOP)
- % KOP, Keyoption. 1--uniform space pattern; 2,3--cosine space pattern
- % N--Total grid point; M--Order
- function [C,X]=DQC(N,M,KOP)
- %1. Generation of grid point coordinates in X-axis
- for i=1:N
- X1(i)=(i-1)/(N-1);
- X2(i)=0.5*(1-cos(pi*(i-1)/(N-1)));
- end
- X3(1)=0;X3(2)=0.00000001;%0.00001;
- for i=3:N
- X3(i)=0.5*(1-cos(pi*(i-1)/(N-1)));
- end
- X4(1)=0;X4(2)=0.00001;X4(N-1)=0.99999;X4(N)=1.0;
- for i=3:N-2
- X4(i)=0.5*(1-cos(pi*(i-2)/(N-3)));
- end
- X5(1)=0;X5(2)=0.00001;X5(3)=0.001;
- X5(N-2)=0.999;X5(N-1)=0.99999;X5(N)=1.0;
- for i=4:N-3
- X5(i)=0.5*(1-cos(pi*(i-3)/(N-5)));
- end
- if KOP==1
- X=X1;
- elseif KOP==2
- X=X2;
- elseif KOP==3
- X=X3;
- elseif KOP==4
- X=X4;
- else
- X=X5;
- end
- %2. Producing the weighting coefficients
- %2.1 Producing the first-order weighting coefficients
- C=zeros(N,N,M);
- for i=1:N
- Q(i)=1.0;
- for j=1:N
- if i~=j
- Q(i)=Q(i)*(X(i)-X(j));
- end
- end
- end
- for i=1:N
- for j=1:N
- if i~=j
- C(i,j,1)=Q(i)/(X(i)-X(j))/Q(j);
- C(i,i,1)=C(i,i,1)-C(i,j,1);
- end
- end
- end
- %2.2 Producing the High-order weighting coefficients
- if M>1
- for k=2:M
- for i=1:N
- for j=1:N
- if i~=j
- C(i,j,k)=k*(C(i,i,k-1)*C(i,j,1)-C(i,j,k-1)/(X(i)-X(j)));
- C(i,i,k)=C(i,i,k)-C(i,j,k);
- end
- end
- end
- end
- end
复制代码
来自:http://hi.baidu.com/hjxiang/blog ... 43902fb93820fb.html |