Matlab 求和中的分支函数

Question

我想计算初始条件 U(x,0)=1, |x| 的单缝衍射的振幅和相位一个，一个=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 和 pre-allocate 矩阵 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

Matlab 求和中的分支函数

Branch function in a sum in Matlab

matlab