Complex equations: NDSolveValue / FEM,

So, I am trying to solve a simple problem using the FEM method.

The distribution of Voltage over a region,

The region is,

  << NDSolveFEM
bm = ToBoundaryMesh[
"Coordinates" -> {{0., 0.}, {1., 0.}, {12/10, 0.}, {12/10, 1.}, {1.,
1.}, {0., 1.}},
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4,
5}, {5, 6}, {6, 1}, {5, 2}}]}];
em = ToElementMesh[bm,
"RegionMarker" -> {{{0.5, 0.5}, 1}, {{1.2, 0.5}, 2}}];


which look like this,

     Show[{
em["Wireframe"["MeshElementStyle" -> Directive[EdgeForm[Red], Thin],
"MeshElementMarkerStyle" -> Blue]],
bm["Wireframe"["MeshElementStyle" -> Black]]
}, ImageSize -> 300]


I have two equations,

op1 = Laplacian[v1[x, y], {x, y}];
op2= Laplacian[v2[x, y], {x, y}];

bc1=DirichletCondition[v1[x, y] == 4, x == 12/10 && y==1];
bc2=DirichletCondition[v2[x, y] == 0, x == 0];

omg=376.991;
dd=100;
d=25*10^-6;

eef={4.18382,4.22462,9.95342,12.8216,13.9473,14.4967,14.8028,14.9896,15.1113,15.195,15.2551,15.2998,15.3341,15.3611,15.3829,15.4007,15.4155,15.428,15.4387,15.4478,15.4556,15.4624,15.4682,15.4733,15.4777};
ee1 = ListInterpolation[eef, {{0, Length[eef] - 1}}, InterpolationOrder -> 1]


ee1 has dependence on v1[x,y]-v2[x,y], as ee1[Abs[v1[x,y]-v2[x,y]]]

{v1fun,v2fun} = NDSolveValue[{op1==If[ElementMarker == 1, ((I omg ee1[
Abs[v1[x, y] - v2[x, y]]])/(dd d)) (v1[x, y] -
v2[x, y]), 0], op2==If[ElementMarker == 1, (
(I omg ee1[
Abs[v1[x, y] - v2[x, y]]])/(dd d)) (v2[x, y] -
v1[x, y]), 0],bc1,bc2},{v1[x, y],v2{x,y}}, {x, y} ∈ em]


The output contour plot of (Abs[v1[x,y]-v2[x,y]]) should be the distribution of voltage over the entire mesh region with the highest values at the x=0 and x=1.2.

I am getting different errors in versions 11.1 and 13.1. I am trying to solve complex equations in NdSolveValue and FEM.

• I think the main problem is that it is not possible to find the derivative of the Abs of complex values. The derivative is needed for linearization. Maybe you can replace the Abs with a smooth function for which a derivative can be found. Feb 10, 2023 at 8:29
• @user21 Thanks! I will follow the suggestion. However, the answer given below is also very close to solution. That could be an another way as well.
– a019
Feb 10, 2023 at 8:33
• Absolutely, I just wanted to express that I am not convinced that the reason is the large numbers, it's a mathematical problem, I think. Feb 10, 2023 at 8:34
• Btw. I added code such that you will not see the Part message any longer (but you will see the other messages). Feb 10, 2023 at 8:36

Parameter k = omg/dd/d of this model is about 150796, together with function ee1 it gives about 10^7 on the right side op1 and op2. It is too high for nonlinear FEM, why NDSolve failed on this example. But the problem could be solved with linear FEM. Therefore we use the false transition method as follows

<< NDSolveFEM
bm = ToBoundaryMesh[
"Coordinates" -> {{0., 0.}, {1., 0.}, {12/10, 0.}, {12/10,
1.}, {1., 1.}, {0., 1.}},
"BoundaryElements" -> {LineElement[{{1, 2}, {2, 3}, {3, 4}, {4,
5}, {5, 6}, {6, 1}, {5, 2}}]}];
em = ToElementMesh[bm,
"RegionMarker" -> {{{0.5, 0.5}, 1}, {{1.2, 0.5}, 2}}];

op1 = Laplacian[v1[x, y], {x, y}];
op2 = Laplacian[v2[x, y], {x, y}];

bc1 = DirichletCondition[v1[x, y] == 4, x == 12/10 && y == 1];
bc2 = DirichletCondition[v2[x, y] == 0, x == 0];

omg = 376.991;
dd = 100;
d = 25*10^-6;

eef = {4.18382, 4.22462, 9.95342, 12.8216, 13.9473, 14.4967, 14.8028,
14.9896, 15.1113, 15.195, 15.2551, 15.2998, 15.3341, 15.3611,
15.3829, 15.4007, 15.4155, 15.428, 15.4387, 15.4478, 15.4556,
15.4624, 15.4682, 15.4733, 15.4777};
ee1 = ListInterpolation[eef, {{0, Length[eef] - 1}},
InterpolationOrder -> 1];
k = omg/dd/d

V1[0][x_, y_] := 4; V2[0][x_, y_] := 0;

Do[{V1[i], V2[i]} =
NDSolveValue[{op1 ==
If[ElementMarker ==
1, ((I k ee1[Abs[V1[i - 1][x, y] - V2[i - 1][x, y]]])) (v1[x,
y] - v2[x, y]), 0],
op2 == If[
ElementMarker ==
1, ((I  k ee1[Abs[V1[i - 1][x, y] - V2[i - 1][x, y]]])) (v2[x,
y] - v1[x, y]), 0], bc1, bc2}, {v1,
v2}, {x, y} \[Element] em];, {i, 5}]


Visualization on every step shows stable solution at i>2

Table[{DensityPlot[Abs[V1[i][x, y]], {x, y} \[Element] em,
ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic, PlotLabel -> "v1"],
DensityPlot[Abs[V2[i][x, y]], {x, y} \[Element] em,
ColorFunction -> "Rainbow", PlotRange -> All,
PlotLegends -> Automatic, PlotLabel -> "v2"]}, {i, 5}]


Update 1. With bc2 = DirichletCondition[v2[x, y] == 0, x == 0&& y==0] we have pictures

• Thanks!! In the problem, V1=4 is at the right up corner of the mesh and V2=0 is at the left down corner of the mesh. However, in density plots, both V1 and V2 at the right up corner. V2 should be at the opposite side of the V1.
– a019
Feb 10, 2023 at 8:28
• @a019 Do you mean that your bc2=DirichletCondition[v2[x, y] == 0, x == 0] is wrong and it should be bc2=DirichletCondition[v2[x, y] == 0, x == 0&&y==0]? Feb 10, 2023 at 10:13
• Yes, I thought x==0 will be enough, however it should be bc2=DirichletCondition[v2[x, y] == 0, x == 0&&y==0]
– a019
Feb 10, 2023 at 10:16
• @a019 See update to my answer. Feb 10, 2023 at 10:20