solvepde fails on complex 'f' coefficient

I am trying to solve a pde using complex coefficients. However I am running into the following issue. Whenever I run solvepde I get an error regarding the complex coefficients.

Here is a MWE:

model = createpde();
geometryFromEdges(model, @circleg);
applyBoundaryCondition(model, 'Edge', 1:model.Geometry.NumEdges, 'u', 0);
generateMesh(model);
fC = @(r, s) complex(1i*ones(1,numel(r.x)));
specifyCoefficients(model, 'm', 0, 'd', 0, 'a', 0, 'c', 1, 'f', fC);
solvepde(model)

Which gives me:

Error using formGlobalKF2D Coefficient evaluation function, "@(r,s)complex(1i*ones(1,numel(r.x)))", returned a complex matrix but previously returned a real matrix. This function must consistently return the same type of matrix.

Any help would be very welcome! Thank you!

Answer 1

[Ravi Kumar]

Hi Yauhen,

The expression

exp((r.x.^2+r.y.^2)/2).*exp(-1i.*s.time)

evaluates to real 1 when s.time =0, r.x= 0, r.y=0. Solver detects switching from real to complex or vice versa. To ensure the value of fC is consistent, please use the following expression for fC.

fC = @(r, s) complex(exp((r.x.^2+r.y.^2)/2).*exp(-1i.*s.time));

In this case fC = 1+0i when s.time =0, r.x= 0, r.y=0. Thus keeping it consistent throughout the solution.

Regards,

Ravi

[Yauhen]

Finally it works!

Thank you very much.

Yauhen

[Hansu zhang]

Thank you Ravi Kumar

My matlab error "returned a complex matrix but previously returned a real matrix. This function must consistently ret". Your expression for fC works well . (ง •̀_•́)ง

Answer 2

[Alan Weiss]

It seems to me that your problem would be best solved by setting the f coefficient to 1i, a constant rather than a function. However, when I tried that, it did not work in solvepde. It is possible that there is a bug in solvepde; the development team is investigating.

Meanwhile, you might do better to use the legacy workflow using assempde.

Sorry, that is all I have at the moment.

Alan Weiss

MATLAB mathematical toolbox documentation

[Avraham Klein]

Actually, placing the 1i constant in specifyCoefficient does work but gives an incorrect answer. This code:

   specifyCoefficients(model, 'm', 0, 'd', 0, 'a', 0, 'c', 1, 'f', 1i);
   solvepde(model)

gives an answer that is identically zero. It appears that solvepde is ignoring the complex part. Giving the following code

   specifyCoefficients(model, 'm', 0, 'd', 0, 'a', 0, 'c', 1, 'f', 1+1i);
   solvepde(model)

yields a solution to the real Poisson equation.

[Yauhen]

Is there any progress on this issue?

[Alan Weiss]

The development team has confirmed that there is a bug in solvepde for some problems with complex coefficients. The solution will be in a future release.

Meanwhile, while the legacy solvers are more cumbersome to use, they do not suffer from this bug.

Sorry about the bug.

Alan Weiss

MATLAB mathematical toolbox documentation

[Yauhen]

Does not work with complex f coefficient in 2018a release. Is there any solution for parabolic PDE with the complex f and d coefficients?

[Ravi Kumar]

Hi Yauhen,

What is the error you get in R2018a. The reproduction as in the original question runs in R2018a. It could be a different issue. Can you provide a reproduction code and error.

Regards, Ravi

[Ravi Kumar]

I also added a dummy complex d-coefficient to the original example, it runs in R2018a.

model = createpde();
geometryFromEdges(model, @circleg);
applyBoundaryCondition(model, 'Edge', 1:model.Geometry.NumEdges, 'u', 0);
generateMesh(model);
fC = @(r, s) complex(1i*ones(1,numel(r.x)));
dC = @(r, s) complex(1i*ones(1,numel(r.x)));
specifyCoefficients(model, 'm', 0, 'd',dC, 'a', 0, 'c', 1, 'f', fC);
setInitialConditions(model,0)
solvepde(model,[0:0.01:0.1])

ans =

TimeDependentResults with properties:

    NodalSolution: [1153×11 double]
    SolutionTimes: [0 0.0100 0.0200 0.0300 0.0400 0.0500 0.0600 0.0700 0.0800 0.0900 0.1000]
       XGradients: [1153×11 double]
       YGradients: [1153×11 double]
       ZGradients: []
             Mesh: [1×1 FEMesh]

[Yauhen]

Hi Ravi Kumar, Thank you for fast reply, here is an example of the code:

model = createpde();
geometryFromEdges(model, @circleg);
applyBoundaryCondition(model, 'Edge', 1:model.Geometry.NumEdges, 'u', 0);
generateMesh(model,'Hmax',0.1);

fC = @(r, s) exp((r.x.^2+r.y.^2)/2).*exp(-1i.*s.time);
dC = -1i;
specifyCoefficients(model, 'm', 0, 'd',dC, 'a', 0, 'c', 1, 'f', fC);
setInitialConditions(model,0);

nframes = 30;
tlist = linspace(0,0.1,nframes);
result = solvepde(model,tlist)
u = result.NodalSolution;

The errors are:

Error using formGlobalKF2D
Coefficient evaluation function, "@(r,s)exp((r.x.^2+r.y.^2)/2).*exp(-1i.*s.time)", returned a complex matrix
but previously returned a real matrix. This function must consistently return the same type of matrix.

Error in pde.DiscretizedPDEModel/getStationaryFEMatrices (line 80)
                [K, F] = formGlobalKF2D(self.emptyPDEModel, self.p, self.t, self.coefstruct,u,time);

Error in pde.DynamicDiscretizedPDEModel/getDynamicFEMatrices (line 118)
            [K,A,F,Q,G,H,R]  = self.getStationaryFEMatrices(u,time);

Error in pde.DynamicDiscretizedPDEModel/checkSpatialCoefsForUorTDependence (line 52)
            [Mass0,K0,A0,F0,Q0,G0,H0,R0] = self.getDynamicFEMatrices(u0, t0);

Error in pde.DynamicDiscretizedPDEModel (line 38)
            obj = obj.checkSpatialCoefsForUorTDependence(u0,tlist);

Error in pde.EquationModel/solveTimeDependent (line 25)
femodel=pde.DynamicDiscretizedPDEModel(self,p,e,t,coefstruct,u0,tlist,tsecondOrder);

Error in pde.PDEModel/solvepde (line 54)
    [u,dudt] = self.solveTimeDependent(coefstruct, u0, ut0, tlist, ...

Do you have any ideas what is wrong here?

Best regards

[Ravi Kumar]

I posted an answer that fixes the error.