Second order derivatives in the PDE Toolbox

The PDE toolbox is pretty brilliant in that it evaluates the spatial derivitives of your function as well as the function itself in PDE results. Say you have an electric field problem, the PDE toolbox can be used to solve for the electric potential, then combined with evaluateGradient, this can be used to calculate the electric field (See 3:68 in the PDE documentation for a 2D example).

My specific problem used the toolbox for a system of 3 equations in three dimensions. The solution is of the form 3xN_nodes, (ux,uy,uz) leading to gradx being an matrix of the same size with components of the form dux/dx, duy/dx, duz/dx.

It is a useful problem then to find the second order derivatives (d2ux/d2x etc) for a whole host of problems e.g finding the turning points in the electric field.

I wondered then if the PDE toolbox (or any other toolbox which can deal with the results of this form) had any tools similar to evaluateGradient for the second order derivatives. I have searched the documentation but didn't see anything useful.

Any suggestions or direction greatly appreciated. I think perhaps the optimisation toolbox would use these derivatives to find maxima or minima but I wouldn't know where to start (or if different toolboxes can even be interfaced).

[Ravi Kumar]

I can suggest something if have a bit more information on your setup. Can you send me an email, I will connect with with you from there.

[Elena Bellizia]

Hi, I have the same question.

An easy example of my setup is the following

model=createpde();
geometryFromEdges(model,@lshapeg);
mesh=generateMesh(model,'Hmax',0.1);
applyBoundaryCondition(model,'dirichlet','edge',1:model.Geometry.NumEdges,'u',0);
specifyCoefficients(model,'m',0,...
    'd',0,...
    'c',1,...
    'a',0,...
    'f',0.5);
results= solvepde(model); 

%Solution at nodes
u = results.NodalSolution;

%Gradient of the solution at nodes
ux=-results.XGradients;
uy=-results.YGradients;

How do I compute second-order derivatives of u?

Thanks

[J. Alex Lee]

What exactly do you mean by "second derivative"? For a scalar field u the "gradient" is a vector field whose components in the reference XY coordinate system are presumably the scalar fields ux and uy you can obtain from results.XGradients and results.YGradients, which are the indepdenet "derivatives" in those orthogoal directions. So if you want to take the gradient again, you'd get another pair of XY derivatives for each of the first X and Y derivatives. That's probably not what you want.

So likely you are both looking for the Laplacian of the scalar field u, which is the divergence of the gradient. Does the pde toolbox give you access to either of these quantities? Maybe it has a field called "Laplacian" or a method to compute "divergence" of an arbitrary vector quantity on the mesh? (Sorry I don't have pde toolbox, nor does it appear super well documented).

But I guess that a problem is that if the toolbox likes to give you solutions only on the nodes (makes sense for visualization purposes), and the underlying method is a second order finite element method (typical, and would make sense since the class of problems TMW can solve looks limited to elliptic scalar equations), you cannot compute unique FEM-consistent approximations to these quantities on the (shared) nodes between elements, because second derivatives on an element are piecewise constant (only first derivatives are continuous on the nodes). That is, each shared node would have different approximations of the Laplacian from the different elements that it belongs to. So there would need to be an extra step to compute the node values from the elemental values by some kind of interpolation or averaging.

Maybe it should be a feature request if it doesn't exist - it does seem like a useful quantity to have access to.