Matlab中求和中的分支函数

Question

我想在初始条件U(x,0)=1， |x|<a和U(x,0)>a, a=2的情况下计算单缝衍射的振幅和相位。我在Matlab中的代码如下

clear all; clc; format long; %2022-BPM based on Fresnel integral
Nx=1024; L=100; x=(L/Nx)*(-Nx/2:Nx/2-1); dx=x(2)-x(1); 
tic  
a=2;
syms u(x);
u(x)=piecewise(-a<x<a, 1, x>a, 0, x<-a, 0);
fplot(u);
S=100; zmax=10; z=linspace(0.1,zmax,S);
for m=1:S
z1=z(m);
for tr=1:Nx
field(tr,m)=exp(1i*pi/4)/(2*1i*sqrt(pi*z1))*...
sum(u.*exp(1i*(x(tr)-x).^2/(4*z1)))*dx;
end;
end; 
[X,Z]=meshgrid(x,z);
figure(1); surf(X',Z',abs(field).^2); shading interp; view([0 90]); 
axis tight; axis square;
figure(2); surf(X',Z',angle(field)); shading interp; view([0 90]); 
axis tight; axis square;
figure(30); plot(x,abs(field(:,S)),'-b*'); hold on;
figure(31); plot(x,angle(field(:,S)),'-b*'); hold on;
I_onaxis=abs(field(Nx/2+1,:)); % on-axis Intensity of the field
figure(32); plot(z,I_onaxis,'-b*'); hold on;
toc

错误在第16行，这里的和是。

Answer 1

不要混合使用符号变量和数字变量。将u(x)函数替换为u = abs(x) < a，预分配矩阵field以提高速度。

clear, clc, format long; %2022-BPM based on Fresnel integral
Nx = 1024; 
L  = 100; 
x  = (L/Nx)*(-Nx/2:Nx/2-1); 
dx = x(2)-x(1); 
tic  
a = 2;
u = abs(x) < a;
S = 100; zmax = 10; z = linspace(0.1,zmax,S);
field = zeros(Nx,S);
for m = 1:S
z1 = z(m);
for tr = 1:Nx
field(tr,m) = exp(1i*pi/4)/(2*1i*sqrt(pi*z1))*...
sum( u .* exp(1i*(x(tr)-x).^2/(4*z1)) ) * dx;
end
end
[X,Z] = meshgrid(x,z);
figure(1); surf(X',Z',abs(field).^2); shading interp; view([0 90]); 
axis tight; axis square;
figure(2); surf(X',Z',angle(field)); shading interp; view([0 90]); 
axis tight; axis square;
figure(30); plot(x,abs(field(:,S)),'-b*'); hold on;
figure(31); plot(x,angle(field(:,S)),'-b*'); hold on;
I_onaxis=abs(field(Nx/2+1,:)); % on-axis Intensity of the field
figure(32); plot(z,I_onaxis,'-b*'); hold on;
toc