function [x,y,Y]=levydistrintheorytest(t0,lambda,T,a,b);

%  计算理论的二维levy分布的图象   画出来的3维曲面图

%   t0 ：levy walk  趋于零的时间

%  lambda ：幂律值

%   T：时间

%   a: 离散的大步长的最大值

%   b：小步长的最大值

error=0.125/180*pi;  %假定平行于坐标轴的误差角度

syms u

c=-gamma(2-lambda)*t0^(lambda-1)*cos(pi*lambda/2)/lambda;

h1=20;  %大的离散步长

h2=0.1; %小的离散步长

x1=-a:h1:-b;

x2=-(b-h2):h2b-h2);

x3=b:h1:a;

x=[x1,x2,x3];

disp('%%%%%%%%%%%计算的次数为%%%%%%%%%%%%%%%')

N=length(x)

y=x;

Z=zeros(N,N);

for i=1:1   %x坐标

i

for j=1:N  % y坐标

theta=angle(complex(x(i),y(j)));

if theta<error && theta>-error  % x轴正半轴

Z(j,i)=double(int(cos(u.*x(i)/(c*T)^(1/lambda))*exp(-(abs(u))^lambda)./pi+cos(1/u.*x(i)/(c*T)^(1/lambda))*exp(-1/(abs(u))^lambda)/u^2./pi,0,1));

elseif     theta<pi/2+error && theta>pi/2-error  % y轴正半轴

Z(j,i)=double(int(cos(u.*y(j)/(c*T)^(1/lambda))*exp(-(abs(u))^lambda)./pi+cos(1/u.*y(j)/(c*T)^(1/lambda))*exp(-1/(abs(u))^lambda)/u^2./pi,0,1));

elseif   theta>pi-error || theta<-pi+error  %x轴负半轴

Z(j,i)=double(int(cos(u.*x(i)/(c*T)^(1/lambda))*exp(-(abs(u))^lambda)./pi+cos(1/u.*x(i)/(c*T)^(1/lambda))*exp(-1/(abs(u))^lambda)/u^2./pi,0,1));

elseif   theta>-pi/2-error && theta<-pi/2+error   %y轴负半轴

Z(j,i)=double(int(cos(u.*y(j)/(c*T)^(1/lambda))*exp(-(abs(u))^lambda)./pi+cos(1/u.*y(j)/(c*T)^(1/lambda))*exp(-1/(abs(u))^lambda)/u^2./pi,0,1));

else

Z(j,i)=double(int(cos(u.*x(i)/(c*T)^(1/lambda)/sin(theta))*exp(-(abs(u))^lambda)./pi+cos(1/u.*x(i)/(c*T)^(1/lambda)/sin(theta))*exp(-1/(abs(u))^lambda)/u^2./pi,0,1))*(double(int(cos(u.*y(j)/(c*T)^(1/lambda)/cos(theta))*exp(-(abs(u))^lambda)./pi+cos(1/u.*y(j)/(c*T)^(1/lambda)/cos(theta))*exp(-1/(abs(u))^lambda)/u^2./pi,0,1)));

%zy=(double(int(cos(u.*y(j)/(c*T)^(1/lambda)/cos(theta))*exp(-(abs(u))^lambda)./pi+cos(1/u.*y(j)/(c*T)^(1/lambda)/cos(theta))*exp(-1/(abs(u))^lambda)/u^2./pi,0,1)));

%Z(j,i)=zx*zy;

end

end

end

Y=Z;

end

错误提示；

Error using mupadmex

Error in MuPAD command: DOUBLE cannot convert the input expression into a

double array.

If the input expression contains a symbolic variable, use the VPA function

instead.

Error in sym/double (line 710)

Xstr = mupadmex('symobj::double', S.s, 0);

Error in levydistrintheorytest (line 37)

Z(j,i)=double(int(cos(u.*x(i)/(c*T)^(1/lambda)/sin(theta))*exp(-(abs(u))^lambda)./pi+cos(1/u.*x(i)/(c*T)^(1/lambda)/sin(theta))*exp(-1/(abs(u))^lambda)/u^2./pi,0,1))*(double(int(cos(u.*y(