MATLAB: Does hyperbolic-solver really use a nonlinear solver?!

Question

My PDE is time-dependent and has a coefficient c , which is depending on the gradient of the solution u (due to the usage of Green-Lagrange-strain measure).

Following the documentation „a nonlinear solver is available for the nonlinear elliptic PDE, the parabolic and hyperbolic equation solvers also solve nonlinear and time-dependent problems.“ For elliptic problems it is necessary to explicitly use pdenonlin instead of assempde . But according to the documentation „the parabolic and hyperbolic functions call the nonlinear solver automatically “.

I don't use the PDE-Tool-GUI, but implemented my own method of line with the pdetool-command-line functions. I looked up the codes of the hyperbolic / parabolic -solver and it's subroutines, and I can nowhere find any call of a nonlinear-solver in case of solution-dependent coefficients!

The results of my calculations also strongly imply, that the nonlinearity induced by c is not treated correctly, since the results are the same as in the linear case. So if anyone can tell me, if it is possible to solve nonlinear transient PDEs with the pdetool-command-line functions, I would be very happy, since I am at my wit's end for another error source producing the unexpected results.

Answer 1

Hi,

Unfortunately, I believe you have uncovered a bug in the hyperbolic function that occurs when you have a nonlinear equation but no Dirichlet boundary conditions. I believe there is a simple fix. If you replace the line

uFull = u;

in toolbox/pde/pdehypdf.m with

uFull = u(1:nu);

your example will run to completion. The line to replace is line 16 in MATLAB R2013b.

However, the results may not be what you intended. I'm guessing that your problem is just a rigid body rotation of the block. If so, the rotated block is larger than the original block which is obviously not correct for the rigid body case.

Before saying more, I should first say that I have wanted to derive a c-matrix for a Lagrangian formulation of 2D elasticity but have been somewhat discouraged by all the algebra required. So I don't have such formulation but believe it is possible to derive one and suspect it is more complicated than what you have in cCoeff_rotation.

I did a simple test of your cCoeff_rotation function (file attached) where I created a simple block, defined the displacement vector for a 45-degree rotation, and calculated the internal forces using cCoeff_rotation. This vector is not zero as expected.

You probably have already looked at this page:

http://www.mathworks.com/help/pde/ug/nonlinear-equations.html?searchHighlight=nonlinear+equations

but there is a key equation there that is easy to overlook. Under the line "The residual vector ρ(U) can be easily computed as" you will see that the residual vector is calculated from a product of K and U. This is different from the typical nonlinear mechanics notation where K is a tangent stiffness matrix.

Hope this documentation page and the fix will point you in the right direction.

Regards,

Bill