一个求方向导数的问题

lyj

nr=44;
ns=29;
Lx=1.0;
Ly=1.5;
c=340;
f=200;
k=2*pi*f/c;
bo=[];
r=0:nr;
k1=r*pi/Lx;                     %求k(r,x),x坐标
k2=sqrt(k^2-k1.^2);           %求k(r,y),y坐标
bo1=cos(k1*x).*exp(-i*k2*y);
s=0:ns;
k1=s*pi/Ly;                  %k(s,y),y坐标
k2=sqrt(k^2-k1.^2);          %k(s,x),x坐标
bo2=exp(-i*k2*x).*cos(k1*y);
bo=[bo1 bo2];                     %波函数矩阵
bo3=bo.';                        %波函数矩阵的转置
x1=0;
y1=0;
x2=1.0;
y2=0;
syms x y
b=sqrt((x2-x1)^2+(y2-y1)^2);
cosa=(x2-x1)/b;
sina=(y2-y1)/b;
A=[cosa sina];      %方向正弦和余弦
bo_direction=[];
r=0:nr;
k1=r*pi/Lx;
k2=sqrt(k^2-k1.^2);
bo1=cos(k1*x).*exp(-i*k2*y);
box1=diff(bo1,x);
boy1=diff(bo1,y);
bon1=[cosa sina]*[box1;boy1];
s=0:ns;
k3=s*pi/Ly;
k4=sqrt(k^2-k3.^2);
bo2=exp(-i*k4*x).*cos(k3*y);
box2=diff(bo2,x);
boy2=diff(bo2,y);
bon2=[cosa sina]*[box2;boy2];
bo_direction=[bon1 bon2];        %波函数的方向导数矩阵
bo_direction1=bo_direction.';
B=bo3*bo_direction;    %波函数的转置与方向导数矩阵相乘
C=bo_direction1*bo;    %方向导数矩阵与波函数矩阵相乘
x3=1.0;
y3=1.5;
x4=0;
y4=0.5;
B1=int(int(B,x,x1,x2),y,y1,y2);
B2=int(int(B,x,x2,x3),y,y2,y3);
B3=int(int(B,x,x3,x4),y,y3,y4);
B4=B1+B2+B3;
C1=int(int(C,x1,x4),y,y3,y4);
a=1.225;
w=400*pi;
b=i/(a*w);
Aaa=b*(B4-C1)

这个怎么出不来啊？帮忙看看

xjzuo

调试了一下你的程序, 发现运行到 B1=int(int(B,x,x1,x2),y,y1,y2);时,符号积分相当耗时.
所以建议改用数值积分.
另: syms x y  应放在程序的前面.

lyj

syms x y
nr=44;
ns=29;
Lx=1.0;
Ly=1.5;
c=340;
f=200;
k=2*pi*f/c;
bo=[];
r=0:nr;
k1=r*pi/Lx;                     %求k(r,x),x坐标
k2=sqrt(k^2-k1.^2);           %求k(r,y),y坐标
bo1=cos(k1*x).*exp(-i*k2*y);
s=0:ns;
k1=s*pi/Ly;                  %k(s,y),y坐标
k2=sqrt(k^2-k1.^2);          %k(s,x),x坐标
bo2=exp(-i*k2*x).*cos(k1*y);
bo=[bo1 bo2];                     %波函数矩阵
bo3=bo.';                        %波函数矩阵的转置
x1=0;
y1=0;
x2=1.0;
y2=0;
b=sqrt((x2-x1)^2+(y2-y1)^2);
cosa=(x2-x1)/b;
sina=(y2-y1)/b;
A=[cosa sina];      %方向正弦和余弦
bo_direction=[];
r=0:nr;
k1=r*pi/Lx;
k2=sqrt(k^2-k1.^2);
bo1=cos(k1*x).*exp(-i*k2*y);
box1=diff(bo1,x);
boy1=diff(bo1,y);
bon1=[cosa sina]*[box1;boy1];
s=0:ns;
k3=s*pi/Ly;
k4=sqrt(k^2-k3.^2);
bo2=exp(-i*k4*x).*cos(k3*y);
box2=diff(bo2,x);
boy2=diff(bo2,y);
bon2=[cosa sina]*[box2;boy2];
bo_direction=[bon1 bon2];        %波函数的方向导数矩阵
bo_direction1=bo_direction.';
B=bo3*bo_direction;    %波函数的转置与方向导数矩阵相乘
C=bo_direction1*bo;    %方向导数矩阵与波函数矩阵相乘
x3=1.0;
y3=1.5;
x4=0;
y4=0.5;
f1=inline('B','x','y');
f2=inline('C','x','y');
B1=dblquad(f1,x1,x2,y1,y2);
B2=dblquad(f1,x2,x3,y2,y3);
B3=dblquad(f1,x3,x4,y3,y4);
B4=B1+B2+B3;
C1=dblquad(f2,x1,x4,y3,y4);
a=1.225;
w=400*pi;
b=i/(a*w);
Aaa=b*(B4-C1)
??? Error using ==> inlineeval
Error in inline expression ==> B
??? Error using ==> eval
Undefined function or variable 'B'.

Error in ==> inline.subsref at 25
    INLINE_OUT_ = inlineeval(INLINE_INPUTS_, INLINE_OBJ_.inputExpr, INLINE_OBJ_.expr);

Error in ==> dblquad>innerintegral at 84
fcl = intfcn(xmin, y(1), varargin{:}); %evaluate only to get the class below

Error in ==> quad at 62
y = f(x, varargin{:});

Error in ==> inline.feval at 20
    [varargout{1:max(1,nargout)}] = builtin('feval',varargin{:});

Error in ==> dblquad at 64
Q = feval(quadf, @innerintegral, ymin, ymax, tol, trace, intfcn, ...

这怎么回事啊？有点蒙了，高手帮忙看看，谢谢

lyj

大家哪位高手帮帮忙啊