11
$\begingroup$

Bug introduced in 10, fixed in 11.1.

I want to solve the electromagnetic wave equation in frequency domain. A known solution, a plane wave Exp[ I k0 x] is used to set the Dirichlet conditions to all the boundaries. I expect to get the plane wave solution. The following code runs well ,but the result is strange. I don't know why. It seems high-frequency oscillation occurs. Is there some options to eliminate the unusual oscillation? thanks a lot! 

<< NDSolve`FEM` λ = 0.53; k0 = 2 π/λ; R = λ; mesh = ToElementMesh[FullRegion[2], {{0, R}, {0, R}}, "MaxCellMeasure" -> 0.0005]; mesh["Wireframe"] op = Most[Curl[Curl[{u[x, y], v[x, y], 0}, {x, y, z}], {x, y, z}] - k0^2 {u[x, y], v[x, y], 0}] pde = op == {0, 0}; Subscript[Γ, D] = DirichletCondition[{u[x, y] == 0., v[x, y] == Exp[I k0 x]}, True]; {us, vs} = NDSolveValue[{pde, Subscript[Γ, D]}, {u, v}, {x, y} ∈ mesh] DensityPlot[Re[vs[x, y]], {x, y} ∈ mesh, ColorFunction -> "Rainbow", PlotLegends -> Automatic, PlotPoints -> 50, PlotRange -> All]

enter image description here

Improve this question
$\endgroup$
10
  • 3
    $\begingroup$ What's your question? $\endgroup$ Commented Oct 23, 2016 at 16:08
  • 4
    $\begingroup$ I'm voting to close this question as off-topic because the code posted does not produce any errors when evaluated. There does not seem to be any question to answer. $\endgroup$ Commented Oct 23, 2016 at 16:37
  • $\begingroup$ The question is probably: "Why does this output not agree with the known analytical solution $v(x,y) = e^{ik_0x}$ and $u(x,y) = 0$?" $\endgroup$ Commented Oct 24, 2016 at 16:34
  • 1
    $\begingroup$ Thank you, I have made a mistake. For a electromagnetic wave , Div(E)==0 should satisfied. To solve Helmholtz equation is the correct choice. $\endgroup$ Commented Oct 24, 2016 at 23:25
  • 1
    $\begingroup$ It seems to work if you replace Exp[I k0 x] by Cos[k0 x] in the DirickletCondition[...] ?? $\endgroup$ Commented Oct 28, 2016 at 20:28

1 Answer 1

Reset to default
9
$\begingroup$

Update

OP's code produces the correct result since v11.1, I think it's safe to say the "noisy" result is caused by a bug.

Here's a solution based on finite difference method (FDM). The definition of pdetoae (a general purpose function discretizing PDE to algebraic equations with finite difference formula) can be found here. 

λ = 53/100; k0 = 2 π/λ; R = λ; op = Most[Curl[Curl[{u[x, y], v[x, y], 0}, {x, y, z}], {x, y, z}] - k0^2 {u[x, y], v[x, y], 0}]; pde = op == {0, 0}; bc = Function[{x, y}, {u[x, y] == 0, v[x, y] == Exp[I k0 x]}] @@@ {{0, y}, {R, y}, {x, 0}, {x, R}} // Flatten; domain = {0, R}; points = 50; grid = Array[# &, points, domain]; difforder = 4; (*Definition of pdetoae isn't included in this code piece, please find it in the link above. *) ptoa = pdetoae[{u, v}[x, y], {grid, grid}, difforder]; del = Most@Rest@# &; ae = del /@ del@# & /@ ptoa@pde; aebc = MapAt[del, ptoa@bc, List /@ Range@4]; {b, mat} = CoefficientArrays[{ae, aebc} // Flatten, Outer[#[#2, #3] &, {u, v}, grid, grid] // Flatten]; sollst = LinearSolve[-mat, N@b]; solmat = ArrayReshape[sollst, {2, points, points}]; {solu, solv} = ListInterpolation[#, {grid, grid}] & /@ solmat; Outer[Plot3D[#@First[#2][x, y], {x, 0, R}, {y, 0, R}, PlotLabel -> #@Last@#2] &, {Re, Im}, {{solu, "u"}, {solv, "v"}}, 1] // Grid

Mathematica graphics

Remark

  1. Number of grid points i.e. points should be even, or the error of $\text{Im}(v)$ will be quite large, I'm not sure about the reason.

  2. Difference order i.e. difforder should be even and greater than 2, or the error will be large and seems unlikely to be alleviated by increasing points. I guess this might explain why OP's code doesn't work well, because according to this document page, the default "MeshOrder" of ToElementMesh is exactly 2, but since I'm still in v9 and the Wolfram Cloud is too slow, I'd like to stop here.

Improve this answer
$\endgroup$
1
  • $\begingroup$ Thank you! you code has help me to understand the EM wave equations better. Considering the ability to handle arbitrarily shaped regions I will continue to study FEM. $\endgroup$ Commented Oct 25, 2016 at 14:12

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.