MATLAB Examples

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),'.-') ```