# Formulation of $-\mathrm{div}(k\nabla u)=f$ in $\Omega$ for the Finite Element Method

% Example1a
%
%   This script demonstrates the use of the Fem code to
%   solve the BVP
%
%                       u=g on Gamma_1,
%                 k*du/dn=h on Gamma_2.
%
%   Here Omega is the unit square, Gamma_1 consists of
%   the top and left edges of the square and Gamma_2 is
%   the remainder of the boundary (the bottom and right
%   edges).  The coefficient is k(x,y)=1+x^2*y, and
%   f, g, h are chosen so that the exact solution is
%   u(x,y)=exp(2*x)*(x^2+y^2).

%   This routine is part of the MATLAB Fem code that
%   accompanies "Understanding and Implementing the Finite
%   Element Method" by Mark S. Gockenbach (copyright SIAM 2006).

% Define the problem functions (note that u and its partial derivatives
% are used here to compute the boundary values; also, these functions
% are used to check the accuracy of the computed solution):

f=vectorize(inline('-2*x*y*(2*exp(2*x)*(x^2+y^2)+2*exp(2*x)*x)-(1+x^2*y)*(4*exp(2*x)*(x^2+y^2)+8*exp(2*x)*x+2*exp(2*x))-2*x^2*exp(2*x)*y-2*(1+x^2*y)*exp(2*x)','x','y'));


This is the first example of some MatLab code of Mark Gockenbach's. He declares that -div(k*grad u)=f. In the first line of MatLab code he declares

f=vectorize(inline('-2*x*y*(2*exp(2*x)*(x^2+y^2)+2*exp(2*x)*x)-(1+x^2*y)*(4*exp(2*x)*(x^2+y^2)+8*exp(2*x)*x+2*exp(2*x))-2*x^2*exp(2*x)*y-2*(1+x^2*y)*exp(2*x)','x','y'));


So there must be some equivalency to -div(k*grad u)=f of such a polynomial. I have a copy of Understanding The Finite Element Method by Gockenbach's to refer to but really did not yet find it where such a derivation is explained. It seems this polynomial is manufactured for the unit square described in the comment block. How is the above polynomial @ ( f=vectorize(inline(... ) derived?

• You could also use the symbolic abilities of matlab to compute f,g,h from the given solution. There might be speed concerns resulting from different levels of interpreted code. Commented Sep 10, 2023 at 16:35

If you take the functions given for $$u$$ and $$k$$ and apply the differential operator that you have in your differential equation you exactly obtain $$f$$, that is

$$f = - \operatorname{div}(k \operatorname{grad}(u))\, .$$

This a method called the Method of Manufactured Solutions and it's used to verify numerical methods. Here, they are assuming some functional forms for $$u$$ and $$k$$ and finding what $$f$$ would correspond to those functions assumed beforehand. See that $$u$$ should satisfy the boundary conditions.