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

7 vues (au cours des 30 derniers jours)
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
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;
%update the plot animation every 5 steps
if mod(round(t/deltat),5)==0
ylim([-3 3])
title([sprintf('t=%f',t*1e6) '\mus'])
%take a snapshot at the user-specified time t=t1
if (t>t1-deltat)&&(t<t1+deltat)
ylim([-3 3])
title([sprintf('t=%f',t*1e6) '\mus'])
%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);

Réponses (2)

Luca Ferro
Luca Ferro le 23 Fév 2023
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 le 23 Fév 2023
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


En savoir plus sur Marine and Underwater Vehicles dans Help Center et File Exchange




Community Treasure Hunt

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

Start Hunting!

Translated by