在使用dblquad时，小弟有用到匿名函数，参照书上写法，小弟的CODE如下:
X=0.01;
C=(2/X)+1;
d1=zeros(C,C);
d2=zeros(C,C);
d3=zeros(C,C);
e1=zeros(C,C);
e2=zeros(C,C);
e3=zeros(C,C);
f1=zeros(9,C^2);
f2=zeros(9,C^2);
f3=zeros(9,C^2);
KAI=1;
ALPHA=4.4;
SIGMA=0.3;
for k2=1:C
k2;
for k1=1:C
k1;
x1=0;
x2=0;
x3=0;
x4=0;
x5=0;
x6=0;
x7=0;
x8=0;
x9=0;
y1=0;
y2=0;
y3=0;
y4=0;
y5=0;
y6=0;
y7=0;
y8=0;
y9=0;
z1=0;
z2=0;
z3=0;
z4=0;
z5=0;
z6=0;
z7=0;
z8=0;
z9=0;
for I=0:20
I;
A=zeros(9,9);
K1=X*(k1)-8-X;
K1;
K2=X*(k2)-8-X;
K2;
A(1,1)=0.5*((K1-pi)^2+(K2-pi)^2)+KAI+(1-SIGMA*i)*(x1*conj(x1)+x2*conj(x2)+2*x3*conj(x3)+2*x4*conj(x4)+2*x5*conj(x5)+2*x6*conj(x6)+2*x7*conj(x7)+2*x8*conj(x8)+2*x9*conj(x9));
A(1,2)=(-0.25*KAI+(1-SIGMA*i)*(0.5*x2*conj(x5)+0.5*x1*conj(x4)+0.5*x3*conj(x6)+0.5*x4*conj(x7)+0.5*x5*conj(x8)+0.5*x6*conj(x9)+0.5*x4*conj(x5)+0.5*x5*conj(x6)+0.5*x7*conj(x8)+0.5*x8*conj(x9)));
A(1,3)=(1-SIGMA*i)*(0.5*x2*conj(x6)+0.5*x4*conj(x6)+x7*conj(x9));
A(1,4)=(-0.25*KAI+(1-SIGMA*i)*(0.5*x4*conj(x7)+0.5*x2*conj(x7)+0.5*x3*conj(x6)+0.5*x2*conj(x5)+x5*conj(x8)+0.5*x6*conj(x9)));
A(1,5)=(1-SIGMA*i)*(0.5*x5*conj(x9)+0.5*x2*conj(x8)+x4*conj(x8)+0.5*x2*conj(x6));
A(1,6)=(1-SIGMA*i)*(0.5*x2*conj(x9)+0.5*x4*conj(x9));
A(1,7)=(1-SIGMA*i)*(x3*conj(x9)+0.5*x2*conj(x8));
A(1,8)=(1-SIGMA*i)*(0.5*x2*conj(x9));
A(1,9)=0;
A(2,1)=-0.25*KAI+(1-SIGMA*i)*(x3*conj(x2)+x5*conj(x4)+x6*conj(x5)+x8*conj(x7)+x9*conj(x8));
A(2,2)=0.5*((K1-pi)^2+(K2)^2)+KAI+(1-SIGMA*i)*(2*x1*conj(x1)+x2*conj(x2)+2*x3*conj(x3)+2*x4*conj(x4)+2*x5*conj(x5)+2*x6*conj(x6)+2*x7*conj(x7)+2*x8*conj(x8)+2*x9*conj(x9));
A(2,3)=-0.25*KAI+(1-SIGMA*i)*(x1*conj(x2)+x7*conj(x8)+x4*conj(x5)+x5*conj(x6)+x7*conj(x8)+x8*conj(x9));
A(2,4)=(1-SIGMA*i)*(x3*conj(x5)+x5*conj(x7)+x6*conj(x8));
A(2,5)=-0.25*KAI+(1-SIGMA*i)*(x1*conj(x4)+x3*conj(x6)+x4*conj(x7)+x1*conj(x4)+x3*conj(x6)+x4*conj(x7)+x5*conj(x8)+0.5*x6*conj(x9));
A(2,6)=(1-SIGMA*i)*(0.5*x1*conj(x5)+0.5*x4*conj(x8)+0.5*x5*conj(x9));
A(2,7)=(1-SIGMA*i)*(1.5*x3*conj(x8));
A(2,8)=(1-SIGMA*i)*(x1*conj(x7)+x3*conj(x9));
A(2,9)=(1-SIGMA*i)*(x1*conj(x8));
A(3,1)=(1-SIGMA*i)*(0.5*x1*conj(x3)+0.5*x4*conj(x6)+0.5*x7*conj(x9)+0.5*x6*conj(x4));
A(3,2)=(-0.25*KAI+(1-SIGMA*i)*(0.5*x2*conj(x3)+0.5*x4*conj(x5)+0.5*x5*conj(x6)+0.5*x7*conj(x8)+0.5*x8*conj(x9)+0.5*x5*conj(x4)));
A(3,3)=0.5*((K1-pi)^2+(K2+pi)^2)+KAI+(1-SIGMA*i)*(x1*conj(x1)+x2*conj(x2)+x3*conj(x3)+2*x4*conj(x4)+2*x5*conj(x5)+2*x6*conj(x6)+2*x7*conj(x7)+2*x8*conj(x8)+2*x9*conj(x9));
A(3,4)=(1-SIGMA*i)*(0.5*x1*conj(x6)+0.5*x2*conj(x5)+x6*conj(x7));
A(3,5)=(1-SIGMA*i)*(x5*conj(x7)+0.5*x2*conj(x6)+0.5*x2*conj(x4)+x6*conj(x8));
A(3,6)=(-0.25*KAI+(1-SIGMA*i)*(0.5*x6*conj(x9)+x4*conj(x7)+x5*conj(x8)+0.5*x1*conj(x4)));
A(3,7)=(1-SIGMA*i)*(0.5*x1*conj(x9)+0.5*x2*conj(x8));
A(3,8)=(1-SIGMA*i)*(0.5*x2*conj(x9));
A(3,9)=0;
A(4,1)=-0.25*KAI+(1-SIGMA*i)*(x5*conj(x2)+x6*conj(x3)+x7*conj(x4)+x8*conj(x5)+x9*conj(x6));
A(4,2)=(1-SIGMA*i)*(x5*conj(x3)+x7*conj(x5)+x8*conj(x6));
A(4,3)=(1-SIGMA*i)*(x7*conj(x6));
A(4,4)=0.5*((K1)^2+(K2-pi)^2)+KAI+(1-SIGMA*i)*(2*x1*conj(x1)+2*x2*conj(x2)+2*x3*conj(x3)+x4*conj(x4)+2*x5*conj(x5)+2*x6*conj(x6)+2*x7*conj(x7)+2*x8*conj(x8)+2*x9*conj(x9));
A(4,5)=-0.25*KAI+(1-SIGMA*i)*(x1*conj(x2)+x2*conj(x3)+x5*conj(x6)+x7*conj(x8)+x8*conj(x9));
A(4,6)=(1-SIGMA*i)*(x1*conj(x3)+0.5*x7*conj(x9));
A(4,7)=-0.25*KAI+(1-SIGMA*i)*(x1*conj(x4)+x2*conj(x5)+x3*conj(x6)+x5*conj(x8)+x6*conj(x9));
A(4,8)=(1-SIGMA*i)*(x1*conj(x5)+x2*conj(x6)+x5*conj(x9));
A(4,9)=(1-SIGMA*i)*(x1*conj(x6));
A(5,1)=(1-SIGMA*i)*(x6*conj(x2)+x8*conj(x4)+x9*conj(x5));
A(5,2)=-0.25*KAI+(1-SIGMA*i)*(x4*conj(x1)+x6*conj(x3)+x7*conj(x4)+x8*conj(x5)+x9*conj(x6));
A(5,3)=(1-SIGMA*i)*(x4*conj(x2)+x7*conj(x5)+x8*conj(x6));
A(5,4)=-0.25*KAI+(1-SIGMA*i)*(x2*conj(x1)+x3*conj(x2)+x6*conj(x5)+x8*conj(x7)+x9*conj(x8));
A(5,5)=0.5*((K1)^2+(K2)^2)+KAI+(1-SIGMA*i)*(2*x1*conj(x1)+2*x2*conj(x2)+2*x3*conj(x3)+2*x4*conj(x4)+x5*conj(x5)+2*x6*conj(x6)+2*x7*conj(x7)+2*x8*conj(x8)+2*x9*conj(x9));
A(5,6)=-0.25*KAI+(1-SIGMA*i)*(x1*conj(x2)+x2*conj(x3)+x4*conj(x5)+x7*conj(x8)+x8*conj(x9));
A(5,7)=(1-SIGMA*i)*(x2*conj(x4)+x3*conj(x5)+x6*conj(x8));
A(5,8)=-0.25*KAI+(1-SIGMA*i)*(x1*conj(x4)+x2*conj(x5)+x3*conj(x6)+x4*conj(x7)+x6*conj(x9));
A(5,9)=(1-SIGMA*i)*(x1*conj(x5)+x2*conj(x6)+x4*conj(x8));
A(6,1)=(1-SIGMA*i)*(x9*conj(x4));
A(6,2)=(1-SIGMA*i)*(x5*conj(x1)+x8*conj(x4)+x9*conj(x5));
A(6,3)=-0.25*KAI+(1-SIGMA*i)*(x4*conj(x1)+x5*conj(x2));
A(6,4)=(1-SIGMA*i)*(x3*conj(x1)+x9*conj(x7));
A(6,5)=-0.25*KAI+(1-SIGMA*i)*(x2*conj(x1)+x3*conj(x2)+x5*conj(x4)+x8*conj(x7)+x9*conj(x8));
A(6,6)=0.5*((K1)^2+(K2+pi)^2)+KAI+(1-SIGMA*i)*(2*x1*conj(x1)+2*x2*conj(x2)+2*x3*conj(x3)+2*x4*conj(x4)+2*x5*conj(x5)+x6*conj(x6)+2*x7*conj(x7)+2*x8*conj(x8)+2*x9*conj(x9));
A(6,7)=(1-SIGMA*i)*(x3*conj(x4));
A(6,8)=(1-SIGMA*i)*(x2*conj(x4)+x3*conj(x5)+x5*conj(x7));
A(6,9)=-0.25*KAI+(1-SIGMA*i)*(x1*conj(x4)+x2*conj(x5)+x3*conj(x6)+x4*conj(x7)+x5*conj(x8));
A(7,1)=(1-SIGMA*i)*(x8*conj(x2)+x9*conj(x3));
A(7,2)=(1-SIGMA*i)*(x8*conj(x3));
A(7,3)=0;
A(7,4)=-0.25*KAI+(1-SIGMA*i)*(x4*conj(x1)+x5*conj(x2)+x6*conj(x3)+x8*conj(x5)+x9*conj(x6));
A(7,5)=(1-SIGMA*i)*(x4*conj(x2)+x5*conj(x3)+x8*conj(x6));
A(7,6)=(1-SIGMA*i)*(x4*conj(x3));
A(7,7)=(0.5*((K1+pi)^2+(K2-pi)^2)+KAI+(1-SIGMA*i)*(2*x1*conj(x1)+2*x2*conj(x2)+2*x3*conj(x3)+2*x4*conj(x4)+2*x5*conj(x5)+2*x6*conj(x6)+x7*conj(x7)+2*x8*conj(x8)+2*x9*conj(x9)));
A(7,8)=-0.25*KAI+(1-SIGMA*i)*(x1*conj(x2)+x2*conj(x3)+x4*conj(x5)+x5*conj(x6)+x8*conj(x9));
A(7,9)=(1-SIGMA*i)*(x1*conj(x3)+x4*conj(x6));
A(8,1)=(1-SIGMA*i)*(x9*conj(x2));
A(8,2)=(1-SIGMA*i)*(x7*conj(x1)+x9*conj(x3));
A(8,3)=(1-SIGMA*i)*(x7*conj(x2));
A(8,4)=(1-SIGMA*i)*(x5*conj(x1)+x6*conj(x2)+x9*conj(x5));
A(8,5)=-0.25*KAI+(1-SIGMA*i)*(x4*conj(x1)+x5*conj(x2)+x6*conj(x3)+x7*conj(x4)+x9*conj(x6));
A(8,6)=(1-SIGMA*i)*(x4*conj(x2)+x5*conj(x3)+x7*conj(x5));
A(8,7)=-0.25*KAI+(1-SIGMA*i)*(x2*conj(x1)+x3*conj(x2)+x5*conj(x4)+x6*conj(x5)+x9*conj(x8));
A(8,8)=0.5*((K1+pi)^2+(K2)^2)+KAI+(1-SIGMA*i)*(2*x1*conj(x1)+2*x2*conj(x2)+2*x3*conj(x3)+2*x4*conj(x4)+2*x5*conj(x5)+2*x6*conj(x6)+2*x7*conj(x7)+x8*conj(x8)+2*x9*conj(x9));
A(8,9)=-0.25*KAI+(1-SIGMA)*(x1*conj(x2)+x2*conj(x3)+x4*conj(x5)+x5*conj(x6)+x7*conj(x8));
A(9,1)=0;
A(9,2)=(1-SIGMA*i)*(x8*conj(x1));
A(9,3)=(1-SIGMA*i)*(x7*conj(x1)+x8*conj(x2));
A(9,4)=(1-SIGMA*i)*(x6*conj(x1));
A(9,5)=(1-SIGMA*i)*(x5*conj(x1)+x6*conj(x2)+x8*conj(x4));
A(9,6)=-0.25*KAI+(1-SIGMA*1i)*(x4*conj(x1)+x5*conj(x2)+x6*conj(x1)+x7*conj(x4)+x8*conj(x5));
A(9,7)=(1-SIGMA*1i)*(x3*conj(x1)+x6*conj(x4));
A(9,8)=-0.25*KAI+(1-SIGMA*1i)*(x2*conj(x1)+x3*conj(x2)+x5*conj(x4)+x6*conj(x5));
A(9,9)=(0.5*((K1+pi)^2+(K2+pi)^2)+KAI+(1-SIGMA*1i)*(2*x1*conj(x1)+2*x2*conj(x2)+2*x3*conj(x3)+2*x4*conj(x4)+2*x5*conj(x5)+2*x6*conj(x6)+2*x7*conj(x7)+2*x8*conj(x8)+x9*conj(x9)));
[ve, ei]=eig(A);
Sum1=0;
Sum2=0;
Sum3=0;
Sum4=0;
Sum5=0;
Sum6=0;
Sum7=0;
Sum8=0;
Sum9=0;
for v=1:1:9;
for u=1:1:9;
Sum1=Sum1+A(v,u)* ve(u,1) *conj(A(v,u)* ve(u,1));
Sum2=Sum2+A(v,u)* ve(u,2) *conj(A(v,u)* ve(u,2));
Sum3=Sum3+A(v,u)* ve(u,3) *conj(A(v,u)* ve(u,3));
Sum4=Sum4+A(v,u)* ve(u,4) *conj(A(v,u)* ve(u,4));
Sum5=Sum5+A(v,u)* ve(u,5) *conj(A(v,u)* ve(u,5));
Sum6=Sum6+A(v,u)* ve(u,6) *conj(A(v,u)* ve(u,6));
Sum7=Sum7+A(v,u)* ve(u,7) *conj(A(v,u)* ve(u,7));
Sum8=Sum8+A(v,u)* ve(u,8) *conj(A(v,u)* ve(u,8));
Sum9=Sum9+A(v,u)* ve(u,9) *conj(A(v,u)* ve(u,9));
end
end
S=[Sum1 Sum2 Sum3 Sum4 Sum5 Sum6 Sum7 Sum8 Sum9];
sorted=sort(S);
m1=sorted(1);
m2=sorted(2);
m3=sorted(3);
index1=find(S<=m1);
index2=find((S>m1)&(S<=m2));
index3=find((S>m2)&(S<=m3));
ei(:,index1);
alpha1=ei(1,index1);
alpha2=ei(2,index1);
alpha3=ei(3,index1);
alpha4=ei(4,index1);
alpha5=ei(5,index1);
alpha6=ei(6,index1);
alpha7=ei(7,index1);
alpha8=ei(8,index1);
alpha9=ei(9,index1);
alpha=alpha1+alpha2+alpha3+alpha4+alpha5+alpha6+alpha7+alpha8+alpha9;
s1=real(alpha);
t1=imag(alpha);
ve(:,index1);
x1=ve(1,index1);
x2=ve(2,index1);
x3=ve(3,index1);
x4=ve(4,index1);
x5=ve(5,index1);
x6=ve(6,index1);
x7=ve(7,index1);
x8=ve(8,index1);
x9=ve(9,index1);
ei(:,index2);
beta1=ei(1,index2);
beta2=ei(2,index2);
beta3=ei(3,index2);
beta4=ei(4,index2);
beta5=ei(5,index2);
beta6=ei(6,index2);
beta7=ei(7,index2);
beta8=ei(8,index2);
beta9=ei(9,index2);
beta=beta1+beta2+beta3+beta4+beta5+beta6+beta7+beta8+beta9;
s2=real(beta);
t2=imag(beta);
ve(:,index2);
y1=ve(1,index2);
y2=ve(2,index2);
y3=ve(3,index2);
y4=ve(4,index2);
y5=ve(5,index2);
y6=ve(6,index2);
y7=ve(7,index2);
y8=ve(8,index2);
y9=ve(9,index2);
ei(:,index3);
kersee1=ei(1,index3);
kersee2=ei(2,index3);
kersee3=ei(3,index3);
kersee4=ei(4,index3);
kersee5=ei(5,index3);
kersee6=ei(6,index3);
kersee7=ei(7,index3);
kersee8=ei(8,index3);
kersee9=ei(9,index3);
kersee=kersee1+kersee2+kersee3+kersee4+kersee5+kersee6+kersee7+kersee8+kersee9;
s3=real(kersee);
t3=imag(kersee);
ve(:,index3);
z1=ve(1,index3);
z2=ve(2,index3);
z3=ve(3,index3);
z4=ve(4,index3);
z5=ve(5,index3);
z6=ve(6,index3);
z7=ve(7,index3);
z8=ve(8,index3);
z9=ve(9,index3);
end
d1(k1,k2)=s1;
e1(k1,k2)=t1;
这边才是用dblquad:
F1=@(x,y)(ALPHA*(x1*exp(i*((K1-pi)*x+(K2-pi)*y))+x2*exp(i*((K1-pi)*x+K2*y))+x3*exp(i*((K1-pi)*x+(K2+pi)*y))+x4*exp(i*((K1)*x+(K2-pi)*y))+x5*exp(i*((K1)*x+(K2)*y))+x6*exp(i*((K1)*x+(K2+pi)*y))+x7*exp(i*((K1+pi)*x+(K2-pi)*y))+x8*exp(i*((K1+pi)*x+(K2)*y))+x9*exp(i*((K1+pi)*x+(K2+pi)*y)*(conj(x1*exp(i*((K1-pi)*x+(K2-pi)*y))+x2*exp(i*((K1-pi)*x+K2*y))+x3*exp(i*((K1-pi)*x+(K2+pi)*y))+x4*exp(i*((K1)*x+(K2-pi)*y))+x5*exp(i*((K1)*x+(K2)*y))+x6*exp(i*((K1)*x+(K2+pi)*y))+x7*exp(i*((K1+pi)*x+(K2-pi)*y))+x8*exp(i*((K1
+pi)*x+(K2)*y))+x9*exp(i*((K1+pi)*x+(K2+pi)*y)).*(-1<=x
& x<=1 & -1<=y & y<=1))))));
F2=@(x,y)(SIGMA*(x1*exp(i*((K1-pi)*x+(K2-pi)*y))+x2*exp(i*((K1-pi)*x+K2*y))+x3*exp(i*((K1-pi)*x+(K2+pi)*y))+x4*exp(i*((K1)*x+(K2-pi)*y))+x5*exp(i*((K1)*x+(K2)*y))+x6*exp(i*((K1)*x+(K2+pi)*y))+x7*exp(i*((K1+pi)*x+(K2-pi)*y))+x8*exp(i*((K1+pi)*x+(K2)*y))+x9*exp(i*((K1+pi)*x+(K2+pi)*y)*(conj(x1*exp(i*((K1-pi)*x+(K2-pi)*y))+x2*exp(i*((K1-pi)*x+K2*y))+x3*exp(i*((K1-pi)*x+(K2+pi)*y))+x4*exp(i*((K1)*x+(K2-pi)*y))+x5*exp(i*((K1)*x+(K2)*y))+x6*exp(i*((K1)*x+(K2+pi)*y))+x7*exp(i*((K1+pi)*x+(K2-pi)*y))+x8*exp(i*((K1
+pi)*x+(K2)*y))+x9*exp(i*((K1+pi)*x+(K2+pi)*y)).*(-1<=x
& x<=1 & -1<=y & y<=1))))));
A1=dblquad(F1,-1,1,-1,1);
A2=dblquad(F2,-1,1,-1,1);
B1=A1/A2;
Ve1=sqrt(B1)*ve(:,index1);
f1(:,k1*k2)=Ve1;
而bug如下:
>> Untitled
Error using *
Inner matrix dimensions must agree.
Error in
@(x,y)(ALPHA*(x1*exp(i*((K1-pi)*x+(K2-pi)*y))+x2*exp(i*((K1-pi)*x+K2*y))+x3*exp(i*((K1-pi)*x+(K2+pi)*y))+x4*exp(i*((K1)*x+(K2-pi)*y))+x5*exp(i*((K1)*x+(K2)*y))+x6*exp(i*((K1)*x+(K2+pi)*y))+x7*exp(i*((K1+pi)*x+(K2-pi)*y))+x8*exp(i*((K1+pi)*x+(K2)*y))+x9*exp(i*((K1+pi)*x+(K2+pi)*y)*(conj(x1*exp(i*((K1-pi)*x+(K2-pi)*y))+x2*exp(i*((K1-pi)*x+K2*y))+x3*exp(i*((K1-pi)*x+(K2+pi)*y))+x4*exp(i*((K1)*x+(K2-pi)*y))+x5*exp(i*((K1)*x+(K2)*y))+x6*exp(i*((K1)*x+(K2+pi)*y))+x7*exp(i*((K1+pi)*x+(K2-pi)*y))+x8*exp(i*((K1+pi
)*x+(K2)*y))+x9*exp(i*((K1+pi)*x+(K2+pi)*y)))))))
Error in quad (line 72)
y = f(x, varargin{:});
Error in dblquad>innerintegral (line 82)
Q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), varargin{:});
Error in quad (line 72)
y = f(x, varargin{:});
Error in dblquad (line 58)
Q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ...
Error in Untitled (line 243)
A1=dblquad(F1,-1,1,-1,1);
请板上高手帮小弟解惑，谢谢:)