8 add one more method for solving the system. edited Feb 20 at 10:01 xzczd 29.5k66 gold badges8383 silver badges273273 bronze badges points = 50; grid = Array[# &, points, domain]; difforder = 2; (*Definition of pdetoae isn't included in this code piece, please find it in the link above.*) ptoa = pdetoae[{Ε["x"], Ε["y"]}[x, y], {grid, grid}, difforder]; del = Most@Rest@# &; ae = del /@ del@# & /@ ptoa@SimplifyPWToUnitStep@pdeS; aebc = MapAt[del, ptoa@bc, List /@ Range@4]; {b, mat} = CoefficientArrays[{ae, aebc} // Flatten, Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten]; sollst = LinearSolve[-mat, N@b]; (* Alternatively, if you're confused about del: *) (* fullae = ptoa@SimplifyPWToUnitStep@pdeS; fullaebc = ptoa@bc; {b, mat} = CoefficientArrays[{fullae, fullaebc} // Flatten, Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten]; sollst = LeastSquares[-mat, N@b, Method -> Direct]; // AbsoluteTiming *) solmat = ArrayReshape[sollst, {2, points, points}]; {solΕx, solΕy} = ListInterpolation[#, {grid, grid}] & /@ solmat; GraphicsGrid[ With[{domain = Sequence @@ (domain + {thickness, -thickness})}, Outer[DensityPlot[#1[First[#2][x, y]], Evaluate[{x, domain}], Evaluate[{y, domain}], PlotRange -> All, PlotLabel -> #1[Last[#2]], ColorFunction -> "AvocadoColors", PlotPoints -> 50] &, {Re, Im}, {{solΕx, "\!$$\*SubscriptBox[\(Ε$$, $$x$$]\)"}, {solΕy, "\!$$\*SubscriptBox[\(Ε$$, $$y$$]\)"}}, 1]], ImageSize -> Large]  points = 50; grid = Array[# &, points, domain]; difforder = 2; (*Definition of pdetoae isn't included in this code piece, please find it in the link above.*) ptoa = pdetoae[{Ε["x"], Ε["y"]}[x, y], {grid, grid}, difforder]; del = Most@Rest@# &; ae = del /@ del@# & /@ ptoa@SimplifyPWToUnitStep@pdeS; aebc = MapAt[del, ptoa@bc, List /@ Range@4]; {b, mat} = CoefficientArrays[{ae, aebc} // Flatten, Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten]; sollst = LinearSolve[-mat, N@b]; solmat = ArrayReshape[sollst, {2, points, points}]; {solΕx, solΕy} = ListInterpolation[#, {grid, grid}] & /@ solmat; GraphicsGrid[ With[{domain = Sequence @@ (domain + {thickness, -thickness})}, Outer[DensityPlot[#1[First[#2][x, y]], Evaluate[{x, domain}], Evaluate[{y, domain}], PlotRange -> All, PlotLabel -> #1[Last[#2]], ColorFunction -> "AvocadoColors", PlotPoints -> 50] &, {Re, Im}, {{solΕx, "\!$$\*SubscriptBox[\(Ε$$, $$x$$]\)"}, {solΕy, "\!$$\*SubscriptBox[\(Ε$$, $$y$$]\)"}}, 1]], ImageSize -> Large]  points = 50; grid = Array[# &, points, domain]; difforder = 2; (*Definition of pdetoae isn't included in this code piece, please find it in the link above.*) ptoa = pdetoae[{Ε["x"], Ε["y"]}[x, y], {grid, grid}, difforder]; del = Most@Rest@# &; ae = del /@ del@# & /@ ptoa@SimplifyPWToUnitStep@pdeS; aebc = MapAt[del, ptoa@bc, List /@ Range@4]; {b, mat} = CoefficientArrays[{ae, aebc} // Flatten, Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten]; sollst = LinearSolve[-mat, N@b]; (* Alternatively, if you're confused about del: *) (* fullae = ptoa@SimplifyPWToUnitStep@pdeS; fullaebc = ptoa@bc; {b, mat} = CoefficientArrays[{fullae, fullaebc} // Flatten, Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten]; sollst = LeastSquares[-mat, N@b, Method -> Direct]; // AbsoluteTiming *) solmat = ArrayReshape[sollst, {2, points, points}]; {solΕx, solΕy} = ListInterpolation[#, {grid, grid}] & /@ solmat; GraphicsGrid[ With[{domain = Sequence @@ (domain + {thickness, -thickness})}, Outer[DensityPlot[#1[First[#2][x, y]], Evaluate[{x, domain}], Evaluate[{y, domain}], PlotRange -> All, PlotLabel -> #1[Last[#2]], ColorFunction -> "AvocadoColors", PlotPoints -> 50] &, {Re, Im}, {{solΕx, "\!$$\*SubscriptBox[\(Ε$$, $$x$$]\)"}, {solΕy, "\!$$\*SubscriptBox[\(Ε$$, $$y$$]\)"}}, 1]], ImageSize -> Large]  7 added 29 characters in body edited Feb 21 '18 at 10:38 xzczd 29.5k66 gold badges8383 silver badges273273 bronze badges Last step is to solve the equation set with finite difference method (FDM). I'll use pdetoae for discretizing:points = 50; grid = Array[# &, points, domain]; difforder = 2; (*Definition of pdetoae isn't included in this code piece, please find it in the link above.*) ptoa = pdetoae[{Ε["x"], Ε["y"]}[x, y], {grid, grid}, difforder]; del = Most@Rest@# &; ae = del /@ del@# & /@ ptoa@SimplifyPWToUnitStep@pdeS; aebc = MapAt[del, ptoa@bc, List /@ Range@4]; {b, mat} = CoefficientArrays[{ae, aebc} // Flatten, Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten]; sollst = LinearSolve[-mat, N@b]; solmat = ArrayReshape[sollst, {2, points, points}]; {solΕx, solΕy} = ListInterpolation[#, {grid, grid}] & /@ solmat; GraphicsGrid[ With[{domain = Sequence @@ (domain + {thickness, -thickness})}, Outer[DensityPlot[#1[First[#2][x, y]], Evaluate[{x, domain}], Evaluate[{y, domain}], PlotRange -> All, PlotLabel -> #1[Last[#2]], ColorFunction -> "AvocadoColors", PlotPoints -> 50] &, {Re, Im}, {{solΕx, "\!$$\*SubscriptBox[\(Ε$$, $$x$$]\)"}, {solΕy, "\!$$\*SubscriptBox[\(Ε$$, $$y$$]\)"}}, 1]], ImageSize -> Large] Remark One may think weWe can solve {pdeS, bc} with NDSolve directly if VersionNumber$VersionNumber >= 1011.1, but sadly it seems not to be true: As one can seeIf you're still in or before v9 (where "FiniteElement" isn't introduced yet), or between $$E_y$$ is full of noise. I believe it's related tov10.0 and v11.0 (where thisa bug pendingisn't fixed yet), finite difference method (FDM) can be used for solving the problem. I'll use pdetoae for discretizing:points = 50; grid = Array[# &, points, domain]; difforder = 2; (*Definition of pdetoae isn't included in this code piece, please find it in the link above.*) ptoa = pdetoae[{Ε["x"], Ε["y"]}[x, y], {grid, grid}, difforder]; del = Most@Rest@# &; ae = del /@ del@# & /@ ptoa@SimplifyPWToUnitStep@pdeS; aebc = MapAt[del, ptoa@bc, List /@ Range@4]; {b, mat} = CoefficientArrays[{ae, aebc} // Flatten, Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten]; sollst = LinearSolve[-mat, N@b]; solmat = ArrayReshape[sollst, {2, points, points}]; {solΕx, solΕy} = ListInterpolation[#, {grid, grid}] & /@ solmat; GraphicsGrid[ With[{domain = Sequence @@ (domain + {thickness, -thickness})}, Outer[DensityPlot[#1[First[#2][x, y]], Evaluate[{x, domain}], Evaluate[{y, domain}], PlotRange -> All, PlotLabel -> #1[Last[#2]], ColorFunction -> "AvocadoColors", PlotPoints -> 50] &, {Re, Im}, {{solΕx, "\!$$\*SubscriptBox[\(Ε$$, $$x$$]\)"}, {solΕy, "\!$$\*SubscriptBox[\(Ε$$, $$y$$]\)"}}, 1]], ImageSize -> Large]  Last step is to solve the equation set with finite difference method (FDM). I'll use pdetoae for discretizing:points = 50; grid = Array[# &, points, domain]; difforder = 2; (*Definition of pdetoae isn't included in this code piece, please find it in the link above.*) ptoa = pdetoae[{Ε["x"], Ε["y"]}[x, y], {grid, grid}, difforder]; del = Most@Rest@# &; ae = del /@ del@# & /@ ptoa@SimplifyPWToUnitStep@pdeS; aebc = MapAt[del, ptoa@bc, List /@ Range@4]; {b, mat} = CoefficientArrays[{ae, aebc} // Flatten, Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten]; sollst = LinearSolve[-mat, N@b]; solmat = ArrayReshape[sollst, {2, points, points}]; {solΕx, solΕy} = ListInterpolation[#, {grid, grid}] & /@ solmat; GraphicsGrid[ With[{domain = Sequence @@ (domain + {thickness, -thickness})}, Outer[DensityPlot[#1[First[#2][x, y]], Evaluate[{x, domain}], Evaluate[{y, domain}], PlotRange -> All, PlotLabel -> #1[Last[#2]], ColorFunction -> "AvocadoColors", PlotPoints -> 50] &, {Re, Im}, {{solΕx, "\!$$\*SubscriptBox[\(Ε$$, $$x$$]\)"}, {solΕy, "\!$$\*SubscriptBox[\(Ε$$, $$y$$]\)"}}, 1]], ImageSize -> Large] Remark One may think we can solve {pdeS, bc} with NDSolve directly if VersionNumber >= 10, but sadly it seems not to be true: As one can see, $$E_y$$ is full of noise. I believe it's related to this pending problem. Last step is to solve the equation set. We can solve {pdeS, bc} with NDSolve directly if $VersionNumber >= 11.1: If you're still in or before v9 (where "FiniteElement" isn't introduced yet), or between v10.0 and v11.0 (where a bug isn't fixed yet), finite difference method (FDM) can be used for solving the problem. I'll use pdetoae for discretizing:points = 50; grid = Array[# &, points, domain]; difforder = 2; (*Definition of pdetoae isn't included in this code piece, please find it in the link above.*) ptoa = pdetoae[{Ε["x"], Ε["y"]}[x, y], {grid, grid}, difforder]; del = Most@Rest@# &; ae = del /@ del@# & /@ ptoa@SimplifyPWToUnitStep@pdeS; aebc = MapAt[del, ptoa@bc, List /@ Range@4]; {b, mat} = CoefficientArrays[{ae, aebc} // Flatten, Outer[#[#2, #3] &, {Ε["x"], Ε["y"]}, grid, grid] // Flatten]; sollst = LinearSolve[-mat, N@b]; solmat = ArrayReshape[sollst, {2, points, points}]; {solΕx, solΕy} = ListInterpolation[#, {grid, grid}] & /@ solmat; GraphicsGrid[ With[{domain = Sequence @@ (domain + {thickness, -thickness})}, Outer[DensityPlot[#1[First[#2][x, y]], Evaluate[{x, domain}], Evaluate[{y, domain}], PlotRange -> All, PlotLabel -> #1[Last[#2]], ColorFunction -> "AvocadoColors", PlotPoints -> 50] &, {Re, Im}, {{solΕx, "\!$$\*SubscriptBox[\(Ε$$, $$x$$]\)"}, {solΕy, "\!$$\*SubscriptBox[\(Ε$$, $$y$$]\)"}}, 1]], ImageSize -> Large] ` 6 edited body edited Sep 20 '17 at 8:44 xzczd 29.5k66 gold badges8383 silver badges273273 bronze badges with Mathematica to make the code more instructive and elegant. The key point is to implementimplementing the generalized curl $$\nabla _s\times$$. There're many possible solutions, for example: with Mathematica to make the code more instructive and elegant. The key point is to implement the generalized curl $$\nabla _s\times$$. There're many possible solutions, for example: with Mathematica to make the code more instructive and elegant. The key point is implementing the generalized curl $$\nabla _s\times$$. There're many possible solutions, for example: 5 edited body edited Apr 30 '17 at 10:56 xzczd 29.5k66 gold badges8383 silver badges273273 bronze badges 4 replaced http://mathematica.stackexchange.com/ with https://mathematica.stackexchange.com/ edited Apr 13 '17 at 12:56 3 added 27 characters in body edited Jan 9 '17 at 3:37 xzczd 29.5k66 gold badges8383 silver badges273273 bronze badges 2 Improve expression a bit. edited Dec 31 '16 at 3:28 xzczd 29.5k66 gold badges8383 silver badges273273 bronze badges 1 answered Dec 30 '16 at 16:44 xzczd 29.5k66 gold badges8383 silver badges273273 bronze badges