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

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

Answer 1

[Luca Ferro]

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.

[Beau Matthews]

Thank you for your answer.

Answer 2

[Walter Roberson]

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

[Beau Matthews]

Thank you for your answer.