# Solve PDEs without transforming to variational formulation **URL:** \ **Category:** variational formulation\ **Created:** [April 29, 2022, 8:22am UTC](https://fenicsproject.discourse.group/t/solve-pdes-without-transforming-to-variational-formulation/8182 "2022-04-29T08:22:56Z")\ **Posts on this page:** 7\ **Page:** 1
**Author:** ![renol-kth](https://avatars.discourse-cdn.com/v4/letter/r/ed655f/32.png) [@renol-kth](https://fenicsproject.discourse.group/u/renol-kth)\ **Post date:** [April 29, 2022, 8:22am UTC](https://fenicsproject.discourse.group/t/solve-pdes-without-transforming-to-variational-formulation/8182/1 "2022-04-29T08:22:56Z")
Dear Community members! I started to learn FEniCS using the book by Langtangen and Logg. When solving the Poisson-equation, it is suggested, to reformulate the equation \Delta u = -f to the variational formulation, i.e. multiply by v \in H^1\_0 and integrate by parts to get \int \nabla u \cdot \nabla v \text{dx} = \int f v \text{dx}. Just to play around with the code, I wanted to implement this equation instead, that one obtains by not doing the integration by parts: \int \Delta u \cdot v \text{dx} = -\int f v \text{dx}. In my understanding u \in H^1 satisfying this expression would still be a weak solution. The code looks like this: ```auto # spaces and mesh mesh = UnitSquareMesh(8, 8) V = FunctionSpace(mesh, 'Lagrange', 1) u = TrialFunction(V) v = TestFunction(V) # BC u_D = Expression('1', degree=0) def boundary(x, on_boundary): return on_boundary bc = DirichletBC(V, u_D, boundary) # equations f = Constant('-10') a = dot(div(grad(u)), v)*dx L = f*v*dx u = Function(V) solve (a==L, u, bc) # plot plot(mesh) plot(u) ``` And produces the following error message: ```auto --------------------------------------------------------------------------- ValueError Traceback (most recent call last) in 21 # plot 22 plot(mesh) ---> 23 plot(u) /usr/local/lib/python3.6/dist-packages/dolfin/common/plotting.py in plot(object, *args, **kwargs) 436 # Plot 437 if backend == "matplotlib": --> 438 return _plot_matplotlib(object, mesh, kwargs) 439 elif backend == "x3dom": 440 return _plot_x3dom(object, kwargs) /usr/local/lib/python3.6/dist-packages/dolfin/common/plotting.py in _plot_matplotlib(obj, mesh, kwargs) 302 303 if isinstance(obj, cpp.function.Function): --> 304 return mplot_function(ax, obj, **kwargs) 305 elif isinstance(obj, cpp.function.Expression): 306 return mplot_expression(ax, obj, mesh, **kwargs) /usr/local/lib/python3.6/dist-packages/dolfin/common/plotting.py in mplot_function(ax, f, **kwargs) 150 if mode == "contourf": 151 levels = kwargs.pop("levels", 40) --> 152 return ax.tricontourf(mesh2triang(mesh), C, levels, **kwargs) 153 elif mode == "color": 154 shading = kwargs.pop("shading", "gouraud") /usr/local/lib/python3.6/dist-packages/matplotlib/tri/tricontour.py in tricontourf(ax, *args, **kwargs) 273 def tricontourf(ax, *args, **kwargs): 274 kwargs['filled'] = True --> 275 return TriContourSet(ax, *args, **kwargs) 276 277 /usr/local/lib/python3.6/dist-packages/matplotlib/tri/tricontour.py in __init__ (self, ax, *args, **kwargs) 33 are described in the docstring of `tricontour`. 34 """ ---> 35 ContourSet. __init__ (self, ax, *args, **kwargs) 36 37 def _process_args(self, *args, **kwargs): /usr/local/lib/python3.6/dist-packages/matplotlib/contour.py in __init__ (self, ax, levels, filled, linewidths, linestyles, alpha, origin, extent, cmap, colors, norm, vmin, vmax, extend, antialiased, *args, **kwargs) 853 self._transform = kwargs.pop('transform', None) 854 --> 855 kwargs = self._process_args(*args, **kwargs) 856 self._process_levels() 857 /usr/local/lib/python3.6/dist-packages/matplotlib/tri/tricontour.py in _process_args(self, *args, **kwargs) 45 else: 46 from matplotlib import _tri ---> 47 tri, z = self._contour_args(args, kwargs) 48 C = _tri.TriContourGenerator(tri.get_cpp_triangulation(), z) 49 self._mins = [tri.x.min(), tri.y.min()] /usr/local/lib/python3.6/dist-packages/matplotlib/tri/tricontour.py in _contour_args(self, args, kwargs) 92 'triangulation') 93 if not np.isfinite(z_check).all(): ---> 94 raise ValueError('z array must not contain non-finite values ' 95 'within the triangulation') 96 ValueError: z array must not contain non-finite values within the triangulation ``` So apparently, the solution has some non-finite values. But, when implementing the left bilinear form as \nabla u \cdot \nabla v the code does produce a solution. I dont understand that behaviour, as the two formulations should be equivalent. Is the problem here, that u \notin \mathbf{C}^2 is somehow assumed in the implementation? I am more than thankful for any kind of help. ---
**Author:** ![mscroggs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mscroggs/32/757_2.png) [@mscroggs](https://fenicsproject.discourse.group/u/mscroggs)\ **Post date:** [April 29, 2022, 11:38am UTC](https://fenicsproject.discourse.group/t/solve-pdes-without-transforming-to-variational-formulation/8182/2 "2022-04-29T11:38:49Z")
I’m guessing some (or maybe all) the coefficients of your solution are nan, which is why you get the error you’re seeing. Degree 1 Lagrange basis functions are continuous between cells, but do not have a continuous derivative. (ie u\in C^1, but u\not\in C^2 as you said). This means that \Delta u is not well defined on the edges of the cells, so to term \displaystyle\int \Delta u\cdot v\,\mathrm{d}x cannot be computed correctly. To see why this is so, you think think about the 1D case. In 1D, the basis functions are hat functions. The derivatives (\nabla u) of these are step functions. The second derivatives (\Delta u) are therefore delta functions. These delta functions should affect the integral (but will not in the FEniCS implementation). In order to assemble \displaystyle\int \Delta u\cdot v\,\mathrm{d}x you will need to use C^2 continuous basis functions. I think these are not supported in FEniCS; and they are not (yet) supported in FEniCSx. ---
**Author:** ![renol-kth](https://avatars.discourse-cdn.com/v4/letter/r/ed655f/32.png) [@renol-kth](https://fenicsproject.discourse.group/u/renol-kth)\ **Post date:** [April 29, 2022, 12:53pm UTC](https://fenicsproject.discourse.group/t/solve-pdes-without-transforming-to-variational-formulation/8182/3 "2022-04-29T12:53:56Z")
Hi Mathew, thank you for your quick answer! So FEniCS does not provide C^2 basis functions - that implies, any equation has to be written in a formulation, where strictly u\in C^1. Does that mean for an equation like this \Delta(\Delta u) = -f one needs to apply integration by parts several times? By using a second order function space ```auto V = FunctionSpace(mesh, 'Lagrange', 2) ``` one would get a basis of second order polynomials, like it is described on your website [DefElement: Lagrange](https://defelement.com/elements/lagrange.html). If the basis functions are second order, they are \in C^2 locally on each element. With these second order elements, I still get the same error from the code, though. Is the problem then caused by lacking regularity of the basis functions on the border/edge of each element? ---
**Author:** ![mscroggs](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/mscroggs/32/757_2.png) [@mscroggs](https://fenicsproject.discourse.group/u/mscroggs)\ **Post date:** [April 29, 2022, 1:09pm UTC](https://fenicsproject.discourse.group/t/solve-pdes-without-transforming-to-variational-formulation/8182/4 "2022-04-29T13:09:52Z")
Hello, Higher order Lagrange elements are still not C^2 continuous, as it is only enforced that their value is the same on mesh entities shared between two cells. In order to have C^2 continuity, you need to have DOFs that endure that the derivative is continuous between cells. For example, Hermite elements ([DefElement: Hermite](https://defelement.com/elements/hermite.html)) are defined using point evaluations and point evaluations of the derivative. If you’re interested, the reason we can’t yet support Hermite elements in FEniCSx is due to needed a more general pull back/push forward map, as described in [this R. Kirby paper](https://doi.org/10.5802/smai-jcm.33). ---
**Author:** ![nate](https://yyz2.discourse-cdn.com/free1/user_avatar/fenicsproject.discourse.group/nate/32/17_2.png) [@nate](https://fenicsproject.discourse.group/u/nate)\ **Post date:** [April 29, 2022, 1:51pm UTC](https://fenicsproject.discourse.group/t/solve-pdes-without-transforming-to-variational-formulation/8182/5 "2022-04-29T13:51:18Z")
In addition to @mscroggs’s points, you could consider weakly enforcing continuity of the first derivative of your FE solution as part of the C^2 requirement. There’s a very old [demo](https://fenicsproject.org/olddocs/dolfin/latest/python/demos/biharmonic/demo_biharmonic.py.html) which could easily be updated for dolfinx. There’s also @kamensky’s [tIGAr](https://github.com/david-kamensky/tIGAr) which offers C^k continuity with a spline basis. E.g. the [bihamonic demo](https://github.com/david-kamensky/tIGAr/blob/master/demos/biharmonic/biharmonic.py). ---
**Author:** ![kamensky](https://avatars.discourse-cdn.com/v4/letter/k/e95f7d/32.png) [@kamensky](https://fenicsproject.discourse.group/u/kamensky)\ **Post date:** [April 29, 2022, 2:11pm UTC](https://fenicsproject.discourse.group/t/solve-pdes-without-transforming-to-variational-formulation/8182/6 "2022-04-29T14:11:33Z")
I’ll also mention an old post [here](https://fenicsproject.discourse.group/t/poisson-equation-the-stubborn-way/1645/3), where I answered a similar question in some more detail. ---
**Author:** ![renol-kth](https://avatars.discourse-cdn.com/v4/letter/r/ed655f/32.png) [@renol-kth](https://fenicsproject.discourse.group/u/renol-kth)\ **Post date:** [May 2, 2022, 7:31am UTC](https://fenicsproject.discourse.group/t/solve-pdes-without-transforming-to-variational-formulation/8182/7 "2022-05-02T07:31:40Z")
Thank you very much everybody!