Skip to main content

1D FDTD Method

· One min read
Kirill Vasin

Solving Maxwell's equation for 1D case with two interfaces. Based on the work of John B. Schneider "Understanding the Finite-Difference Time-Domain Method" 2023 we will try to solve 1D Maxwell's equation for electro-magnetic wave propagating through a non-dispersive plate

Download original notebook

Define constants

ζ = 377.0; (*BB[*)(*Ohm*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
size = 400;

ez = ConstantArray[0., {size}];
hy = ConstantArray[0., {size}];

cμ = ConstantArray[ζ,   {size}];
cϵ = ConstantArray[1.,  {size}];

(*BB[*)(*Define a second medium with losses*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)

cϵ[[200;;220]] = (*FB[*)(((1 - 0.001))(*,*)/(*,*)((1 + 0.001)))(*]FB*);
cμ[[200;;220]] = (*FB[*)((ζ)(*,*)/(*,*)(9.0))(*]FB*) / (1 + 0.001);

generator = With[{ie = size}, Compile[{{steps, _Integer}},
    Module[
     {ez = ez, hy = hy},
     Do[

        (*BB[*)(*Calculate magnetic field *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
        Do[
           hy[[j]] += (*FB[*)(((ez[[j + 1]] - ez[[j]]))(*,*)/(*,*)(ζ))(*]FB*);
        , {j, 1, ie-1}];

        (*BB[*)(*Boundary conditions*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
        hy[[ie]] = hy[[ie-1]];
        ez[[ie]] = ez[[ie-1]];
      
        (*BB[*)(*Pulse generator*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
        hy[[49]] -= Exp[- (*FB[*)(((*SpB[*)Power[(n - 30.)(*|*),(*|*)2](*]SpB*))(*,*)/(*,*)(100.))(*]FB*)] / ζ;

        ez[[1]] = ez[[2]];

        (*BB[*)(*Calculate electric field *)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
        Do[
          ez[[j]] = cϵ[[j]] ez[[j]] + cμ[[j]] (hy[[j]] - hy[[j - 1]])
        , {j, 2, ie - 1}];

        (*BB[*)(*Pulse generator*)(*,*)(*"1:eJxTTMoPSmNhYGAo5gcSAUX5ZZkpqSn+BSWZ+XnFaYwgCS4g4Zyfm5uaV+KUXxEMUqxsbm6exgSSBPGCSnNSg9mAjOCSosy8dLBYSFFpKpoKkDkeqYkpEFXBILO1sCgJSczMQVYCAOFrJEU="*)(*]BB*)
        ez[[50]] += Exp[- (*FB[*)(((*SpB[*)Power[(n + 0.5 - (-0.5) - 30.)(*|*),(*|*)2](*]SpB*))(*,*)/(*,*)(100.0))(*]FB*)];

        
        
      , {n, steps}]; ez],
    CompilationOptions -> {"InlineExternalDefinitions" -> True}, 
    "CompilationTarget" -> "C", "RuntimeOptions" -> "Speed"]];

Plot as a function of x-coordinate

electricField = {Range[400], generator[3]} // Transpose;
Graphics[
  Line[electricField // Offload], 
  PlotRange->{{0,400}, {-1,1}}, 
  TransitionDuration->5
]
(*VB[*)(FrontEndRef["c163d4ea-9002-4a25-a9f0-9e37f0a858d1"])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKJxuaGaeYpCbqWhoYGOmaJBqZ6iZaphnoWqYam6cZJFqYWqQYAgCBDhVu"*)(*]VB*)

Try to drag this slider 😉

Manipulate the time

EventHandler[InputRange[1,800,20,1], Function[val, 
  electricField = {Range[400], generator[val]} // Transpose;
]]
(*VB[*)(EventObject[<|"Id" -> "97462fa3-8bc0-472a-919b-54cc650eacf1", "Initial" -> 1, "View" -> "f294fa8d-e0cd-416c-9be1-baffe32749b6"|>])(*,*)(*"1:eJxTTMoPSmNkYGAoZgESHvk5KRCeEJBwK8rPK3HNS3GtSE0uLUlMykkNVgEKpxlZmqQlWqTophokp+iaGJol61ompRrqJiWmpaUaG5mbWCaZAQCVOxZr"*)(*]VB*)

Congratulations! 🥂 Now you know how to solve 1D Maxwell's equation with an arbitary boundary conditions!