Unexpected results obtained for Poisson problem in complex function space (R2017, R2018)

I'm using MATLAB 2017a. I want to solve Poisson problem with PDEToolbox

div grad u = f

with mixed boundary conditions, where both f and boundary conditions can be complex.

My MWE is given below and it is intended to solve the problem on unit square domain with

f = 1 + i

u = 1 + i on top,left and bottom edges of the square and

du//dn = 1+i on the right edge of the square.

Since real and imaginary parts of BCs and f are the same I was expecting both real and imaginary parts of the solution to be equal. Unfortunately, I keep obtaining different plots for real and imaginary parts of the solution -> see figures below (Figure 2 and Figure 3 produced by MWE).

Even simple numerical differentiation reveals that for imaginary part Neumann condition of the right edge is not fulfilled (lines 60-65, Figures 4 and 6 produced by MWE). It's not equal 1 as it should be, but -1.

I discovered, that I need to impose du//dn = 1-i on the right edge (line 30 instead of line 28) to obtain imaginary part of the solution the same as real, but again, in this way, Neumann BC isn't still fulfilled.

Is it a Matlab's bug or maybe is there some math involved here that I'm missing?

function MWE()

g = @UnitSquareg;
c = 1;
a = 0;
f = 'ones(size(x)) + i*ones(size(x))' ;

%% Create a PDE Model with a single dependent variable
numberOfPDE = 1;
pdem = createpde(numberOfPDE);

%% Create the geometry and append to the PDE Model
geometryFromEdges(pdem,g);

%% Boundary Conditions
%%% Plot the geometry and display the edge labels for use in the boundary condition definition.
figure; 
pdegplot(pdem, 'edgeLabels', 'on'); 
axis([-0.1 1.1 -0.1 1.1]);
title('Geometry With Edge Labels Displayed')

% applying boundary conditions
%       - Dirichlet: h*u=r
%       - Neumann: n*c*grad(u)+qu = g
b1 = applyBoundaryCondition(pdem,'dirichlet','Edge',[1 3 4], 'u', 1+i*1);
% bR = applyBoundaryCondition(pdem,'dirichlet', 'Edge',[ 2 ], 'u', -(1+i*1));

bR = applyBoundaryCondition(pdem,'neumann', 'Edge',[ 2 ], 'q', 0, 'g', (1+i*1));
% uncomment the line below to obtain the same real and imaginary part of the solution
% bR = applyBoundaryCondition(pdem,'neumann', 'Edge',[ 2 ], 'q', 0, 'g', (1-i*1));

%% Create Mesh
hMax = 0.05 ;
mesh = generateMesh(pdem, 'Hmax', hMax, 'GeometricOrder', 'linear');

%% Solution 
u = assempde(pdem,c,a,f);

%% Interpolation of the solution
[p,e,t] = meshToPet(pdem.Mesh) ;
pdeInt = pdeInterpolant(p,t,u) ;
uInterp = @(x,y) reshape( pdeInt.evaluate(x, y), size(x,1), size(x,2) ) ;

%% Prepare mesh for plotting
x = linspace(0,1,41) ;
y = linspace(0,1,41) ;
dx = x(2)-x(1) ;
dy = y(2)-y(1) ;
[X,Y] = meshgrid(x,y) ;

figure,
surf(X,Y, real(uInterp(X,Y)))
colormap jet
xlabel('x')  ; ylabel('y')  ; title('real')
figure,
surf(X,Y, imag(uInterp(X,Y)))
colormap jet
xlabel('x')  ; ylabel('y')  ; title('imag')

[FX,FY] = gradient(uInterp(X,Y), dx,dy) ;
figure, surfc(X,Y,real(FX) ); title('Re Flux X')
figure, surfc(X,Y,real(FY) ); title('Re Flux Y')

figure, surfc(X,Y,imag(FX) ); title('Im Flux X')
figure, surfc(X,Y,imag(FY) ); title('Im Flux Y')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function [x,y]=UnitSquareg(bs,s)
%SQUAREG   Gives geometry data for the squareg PDE model

xmin = 0 ; xmax = 1 ;  ymin = 0 ; ymax = 1 ;

nbs=4;

if nargin==0
  x=nbs; % number of boundary segments
  return
end

d=[
  0 0 0 0 % start parameter value
  1 1 1 1 % end parameter value
  0 0 0 0 % left hand region
  1 1 1 1 % right hand region
];

bs1=bs(:)';

if find(bs1<1 | bs1>nbs)
  error(message('pde:squareg:InvalidBs'))
end

if nargin==1
  x=d(:,bs1);
  return
end

x=zeros(size(s));
y=zeros(size(s));
[m,n]=size(bs);
if m==1 && n==1
  bs=bs*ones(size(s)); % expand bs
elseif m~=size(s,1) || n~=size(s,2)
  error(message('pde:squareg:SizeBs'));
end

if ~isempty(s)

% boundary segment 1
ii=find(bs==1);
if length(ii)
x(ii)=interp1([d(1,1),d(2,1)],[xmin xmax],s(ii));
y(ii)=interp1([d(1,1),d(2,1)],[ymax ymax],s(ii));
end

% boundary segment 2
ii=find(bs==2);
if length(ii)
x(ii)=interp1([d(1,2),d(2,2)],[xmax xmax],s(ii));
y(ii)=interp1([d(1,2),d(2,2)],[ymax ymin],s(ii));
end

% boundary segment 3
ii=find(bs==3);
if length(ii)
x(ii)=interp1([d(1,3),d(2,3)],[xmax xmin],s(ii));
y(ii)=interp1([d(1,3),d(2,3)],[ymin ymin],s(ii));
end

% boundary segment 4
ii=find(bs==4);
if length(ii)
x(ii)=interp1([d(1,4),d(2,4)],[xmin xmin],s(ii));
y(ii)=interp1([d(1,4),d(2,4)],[ymin ymax],s(ii));
end

end

[David Goodmanson]

Hi tri,

I don't have the pde toolbox, but it may solve for -del^2 = f, not del^2= f. That issue aside, mathematically the result you are expecting is correct. What happens when you replace 1+i with i everywhere?

[tri stan]

Hi David,

a) solving -del^2 = f instead of del^2= f results in solutions that are inversions of the figures above,

b) If I replace 1+i with i everywhere I get real part of the solution is equal zero, while imaginary part stays the same.

Edit:

I had a chance to run my code on Matlab 2018. I had to modify the code a bit (because generateMesh function from R2018a by default creates geometry with quadratic elements, while assempde still requires them to be linear!), yet now code should run flawlessly both on R2017 and R2018.

Nevertheless, the results are still the same (erroneous?).