www.pudn.com > chearray.rar > chearray.m, change:2017-03-08,size:1887b


clear all; 
syms  x  
N=8; 
SLL=26; 
r0=20;%10^(SLL/20); 
x0=1.15;%(1/2)*((r0+sqrt(r0^2-1))^(1/(N-1))+(r0-sqrt(r0^2-1))^(1/(N-1)));  
if mod(N,2)==0  %创建A矩阵 
    m=N/2; 
    j=1; 
 A = sym(zeros(m,j)); 
for i=1:m 
        cmd = sprintf('sym(''A_%i_%i'')',i,j); 
        A(i,j) = eval(cmd); 
end 
        A; 
  g00=1;g01=x/x0;f00=x/x0; 
     f(1,1)=f00; 
     for k=1:1:N-2 
         g=2*(x/x0).*g01-g00; 
         g00=g01; 
         g01=g; 
         if mod(k,2)==0 
            f(1,k/2+1)= g01; 
         end 
     end 
     E8=f*A;%计算总场强 
     for p=1:2:N-1 
         a(1,(p+1)/2)=diff(E8,x,p); 
         b(1,(p+1)/2)=a(1,(p+1)/2)/factorial(p); 
         d(1,(p+1)/2)=subs(b(1,(p+1)/2),'x',0); 
     end 
  for i=1:1:m 
         c=d(1,i); 
         for j=1:1:m 
             C(i,j)=diff(c,A(j,1)); 
         end 
  end 
     C 
else 
    m=(N-1)/2; %N为奇数的算法 
    j=1; 
 A = sym(zeros(m+1,j)); 
for i=1:m+1 
        cmd = sprintf('sym(''A_%i_%i'')',i,1); 
        A(i,1) = eval(cmd); 
end 
   A 
   g00=1;g01=x/x0;f00=1; 
     f(1,1)=f00; 
      i=1; 
     for k=1:1:N-2 
         g=2*(x/x0).*g01-g00; 
         g00=g01; 
         g01=g; 
         if mod(k,2)==1 
             i=i+1; 
           f(1,i)=g01 
         end 
     end 
     E8=f*A; 
     d(1,1)=A(1,1); 
     for p=4:2:N+1 
         a(1,p/2)=diff(E8,x,p-2); 
         b(1,p/2)=a(1,p/2)/factorial(p-2); 
         d(1,p/2)=subs(b(1,p/2),'x',0); 
     end 
  for i=1:1:m+1 
         c=d(1,i); 
         for j=1:1:m+1 
             C(i,j)=diff(c,A(j,1)); 
         end 
     end 
end 
   T00=1;T01=x; 
     for p=1:1:N-2 %计算切比雪夫多项式 
         T=2*x.*T01-T00; 
         T00=T01; 
         T01=T; 
    end 
     T1=T01; 
      B=coeffs(T1);  %令场强表达式的系数与切比雪夫多项式系数相等 
             D=B'; 
             A=inv(C)*D 
             P=vpa(A) 
             m=max(P); 
             Ak=m.\P