I am having a hard time getting the matrix multiplication to work due to incompatible matrix sizes.

3 次查看（过去 30 天）
Beau Matthews
Beau Matthews 2023-2-23
评论： Beau Matthews 2023-2-23
I am having trouble with getting the correct matrix multiplication sizes correct at line 35 for the 1-D scalar update equation. Any help would be much appreciated.
% This code implements a one-dimensional scalar wave equation following the
% finite-difference time domain (FDTD) method of Taflove and Hagness
% Physical constants and user-defined values
mu_0 = pi*4e-7; %permeability of free space in H/m
eps_0 = 8.854e-12; %permittivity of free space in F/m
eps_r = 10; %relative permittivity of the material
c = 1/sqrt(mu_0*eps_0); %speed of light in free space in m/s
f = 20*(10^6) ; %frequency of the input wave in Hz
lambda = c/f; %free space wavelength in m
% Set up the space and time dimensions
deltax = 45; %user-defined spatial step in m
xmax = 900; %max value of distance (x) in the simulation space
x = (0:deltax:xmax); %all values of distance (x) in the simulation space
tmax = 2.7*(10^-4); %simulation stops after time t=tmax (in seconds)
S = 2; %user-defined Courant stability factor
deltat = S*deltax/c; %time step in seconds
t=0; %start time in seconds
% Put dielectric material into the simulation space
eps = ones(length(x),1).*eps_0; %initialize permittivity everywhere
s1 = xmax/2; %location index of the material boundary
eps(s1:end) = eps(s1:end)*eps_r;%add dielectric material past boundary
% Initialize electric field
E = zeros(length(x),3); %E=0 everywhere and for all previous time
% Set locations for the virtual electric field probes
x0 = 1; %index of field probe E0
x1 = 451; %index of field probe E1
t1 = 1.35*(10^-4); %time of the snapshot in seconds
%create a blank figure for the FDTD animation
h=figure;
dummy = ((1/(mu_0*eps_0.*eps)).*((deltat/deltax)^2).*((E(3:end,2))-(2*E(2:end-1,2))+(E(1:end-2,2))+(2*E(2:end-1,2))-(E(2:end-1,1))));
% Begin the FDTD update loop
while t<tmax %update until the max time value is reached
%implement the 1-D scalar update equation
E(2:end-1,3) = ((1/(mu_0*eps_0.*eps)).*((deltat/deltax)^2).*((E(3:end,2))-(2*E(2:end-1,2))+(E(1:end-2,2))+(2*E(2:end-1,2))-(E(2:end-1,1))));
%E-field of the incoming wave (turns off after 5 cycles)
if t<=5/f
E(1,3) = sin(2*pi*f*t);
elseif t>5/f
E(1,3) = 0;
end
%update the plot animation every 5 steps
if mod(round(t/deltat),5)==0
figure(h)
plot(x,E(:,3))
xline(x(s1))
ylim([-3 3])
title([sprintf('t=%f',t*1e6) '\mus'])
drawnow
end
%take a snapshot at the user-specified time t=t1
if (t>t1-deltat)&&(t<t1+deltat)
figure
plot(x,E(:,3))
xline(x(s1))
ylim([-3 3])
title([sprintf('t=%f',t*1e6) '\mus'])
end
%store values for the virtual field probes
E0(round(t/deltat) + 1) = E(x0,3);
E1(round(t/deltat) + 1) = E(x1,3);
%move forward one time step
t = t+deltat;
E(:,1) = E(:,2);
E(:,2) = E(:,3);
end

请先登录，再回答此问题。

回答（2 个）

Luca Ferro
Luca Ferro 2023-2-23
I found the issue, i'll guide you through it but unfortunatly i cannot solve it for you since i'm not sure on what your goal is.
E(2:end-1,3) = ((1/(mu_0*eps_0.*eps)).*((deltat/deltax).^2).*((E(3:end,2))-(2*E(2:end-1,2))+(E(1:end-2,2))+(2*E(2:end-1,2))-(E(2:end-1,1))));
in this line the left term has size 19x1 while the right 19x21. Clearly you cannot perform this assignment.
Upon further inspection on the right term, (1/(mu_0*eps_0.*eps) is what causes the unwanted change of dimensions since eps is 21x1 while everything else ((deltat/deltax)^2).*((E(3:end,2))-(2*E(2:end-1,2))+(E(1:end-2,2))+(2*E(2:end-1,2))-(E(2:end-1,1))) is 19x1, meaning that their multiplication will result in a 19x21.
Going up the code eps is initialized as eps = ones(length(x),1).*eps_0 so it is depending on x, which is the root of the problem since it has dimension 1x21 being declared as x=0:deltax:xmax.
What you have to do is change deltax and xmax so that x is of the right dimensions.
Walter Roberson
Walter Roberson 2023-2-23
You have problems with the / operation. / is matrix division with A/B being similar to A * pinv(B)
You need element by element division which is the ./ operation.
As a guideline, I recommend that you do not use the / operator unless it is / by a scalar constant. If you have an actual matrix division A/B rewrite it in terms of (B' \ A')' . Avoiding using / will prevent the kind of problem that you are seeing where 1/expression is not giving you the result you expect

类别

在 Help Center 和 File Exchange 中查找有关 Electromagnetism 的更多信息

标签

产品


版本

R2022b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by