Maxwell's equations

Toby Driscoll, November 10, 2010

(Chebfun example pde/Maxwell.m)

In one dimension in a simple medium, Maxwell's equations can be written as

E_t = c^2 B_x,   B_t = E_x

where E is the z-component of the electric field, B is the y-component of the magnetic field, and c is the speed of light in the medium. A typical boundary condition is to set E=0, which implies a perfect conductor.

Here is how we can pose Maxwell's equations for solving in PDE15S, in a medium with c=1.

pdefun = @(E,B,diff) [ diff(B), diff(E) ];
bc.left = @(E,B,diff) E;
bc.right = bc.left;

Note how the spatial differentiation operator has to be accepted as an input argument in each case. Now we can define an initial condition and solve the problem.

d = [-2,2];
x = chebfun('x',d);
E0 = exp(-16*x.^2);
t = 0:0.05:5;
soln = pde15s(pdefun,t,[E0,-E0],bc,pdeset('plot','on'));

The result is a cell array: soln{1} is a quasimatrix with each column representing E at a requested time, and soln{2} is the same for B. A space-time plot shows an inverted reflection of E off of the right boundary:

[E,B] = deal(soln{:});
waterfall(E,t,'simple')
view(10,72)

In theory, energy is conserved. We can use this principle to check the accuracy of the numerical solution.

energy = sum( E.^2 + B.^2, 1 );  % integrate E^2 + B^2 in space
plot(t,energy-energy(1),'.-')

Clearly, a little dissipation is introduced by the time integration. The amount is compatible with the default tolerance of 1e-6 in PDE15S.

Because the problem is a linear evolution u_t = Au, where u=[E,B], we can use an operator exponential as an alternative for the time evolution. The block operator A is defined as

A = chebop(d);
A.op = @(x,E,B) [ diff(B), diff(E) ];
A.lbc = @(E,B) E;
A.rbc = @(E,B) E;

Now we use expm to solve the linear evolution. We'll save time if we just apply one exponentiation repeatedly.

Phi = expm( (t(2)-t(1))*A );
u = [E0,-E0];
for n = 2:length(t), u=Phi*u; E(:,n) = u(:,1); B(:,n) = u(:,2);
plot(u), drawnow, pause(eps), end

For linear propagation, expm is far more accurate than pde15s, as a check of energy conservation confirms.

energy = sum( E.^2 + B.^2, 1 );  % integrate E^2 + B^2 in space
plot(t,energy-energy(1),'.-')

Published with MATLAB® 7.14