Nonlinear FEM and FindRoot

I'm trying to develop a kind of nonlinear FEM application using mathematica to solve a bvp like the following:

$$\gamma(u') ~u^{iv} + 2 \gamma'(u') u''' u''+ u''^3 = f(x)$$ where $u = \tilde{u}(x)$, $\gamma = \tilde{\gamma}(u')$ and $'$ denotes derivative with respect to the function variable, thus:

$$u' = \frac{\text{d} u(x)}{\text{d}x} \quad \quad \gamma' = \frac{\text{d} \gamma(u')}{\text{d} u'}$$

The next step is getting the weak form:

$$\int_{\Omega} \left( \gamma ~u^{iv} + 2 \gamma' u''' u''+ u''^3 \right) \psi = \int_{\Omega} f(x)\psi$$

where $\Omega$ is the domain and $\psi = \tilde{\psi}(x)$ is the test function. In my code, I'm using a Cubic Hermite Shape Functions, then the approximate function reads as follows:

$$\hat{u} = \sum_i \alpha_i \psi_i^{(a)} + \beta_i \psi_i^{(b)}$$

where $\psi_i^{(a)}$ and $\psi_i^{(b)}$ are piecewise functions, and $\alpha_i$ and $\beta_i$ are coefficients related to the values of function and its first derivative, respectively. Let bear in mind the chain rule:

$$\frac{\text{d} \gamma(u')}{\text{d} x} = \frac{\text{d} \gamma(u')}{\text{d} u'} \frac{\text{d} u'}{\text{d} x} = \gamma' u''$$

After some calculation, the weak form assumes the following form:

$$\int_{\Omega}\gamma~u^{iv} = \left[\left(\gamma u''' + \frac{1}{2}\gamma' u''^2 \right)\psi-\left(\gamma \psi \right)' u'' \right]_{\partial \Omega} + \int_{\Omega} \frac{1}{2} \gamma'' u''^3 \psi + \frac{3}{2} \gamma' u''^2 \psi' + \gamma u'' \psi ''$$

$$\int_{\Omega} 2 \gamma' u''' u'' \psi = \left[\gamma' u''^2 \psi \right]_{\partial \Omega} - \int_{\Omega} \gamma'' \psi u''^3 \psi + \gamma' u''^2 \psi'$$

As starting point, however, I'm using a linear equation, setting $\gamma(u') = 1$ (but at the end, the procedure is the same). Then, the simple equation is:

$$u^{iv} = f(x) \Longrightarrow \int_{\Omega} u^{iv} \psi = \int_{\Omega} f(x) \psi$$ then $$\left[u''' \psi - u'' \psi' \right]_{\partial \Omega} + \int_{\Omega} u'' \psi'' - \int_{\Omega} f(x) \psi = 0$$ If the domain is discretised in n points (n-1 elements), this equation represents a system of 2n (non)linear equations. Let suppose to apply BCs on the function and its first derivative at both edges. The system reads as follows:

$$R(1) = \alpha_1 - u(x_0) \\\\ R(2) = \beta_1 - u'(x_0) \\\\ ... \\\\ R(i) = \int_{\Omega_i} u'' \psi'' - \int_{\Omega_i} f(x) \psi \\\\ ... \\\\ R(2n-1) = \alpha_n - u(x_1) \\\\ R(2n) = \beta_n - u'(x_1) \\\\$$

where the domain is $\Omega = [x_0, x_1]$. Now the goal is solving this system (in general nonlinear). I'm trying to do that using FindRoot. Thus I defined the following functions:

(* integration of forcing term *)
funFF[X_, f_] := Module[
{F, n},
n = Length[X];
F = ConstantArray[0, 2 n];
Do[
jr = 2 j - 1;
F[[jr ;; jr + 3]] =
F[[jr ;; jr + 3]] +
Flatten[NIntegrate[
f[x] NN0[x, X[[j]], X[[j + 1]]], {x, X[[j]], X[[j + 1]]}]],
{j, 1, n - 1}
];
F
];

(* assembly of system *)
SYS1[X_, U_?VectorQ, RHS_, fun_] := Module[
{EQ, x1, x2, f1, p1, f2, p2, F, P, R, n},
n = Length[X];
(* value of function *)
F = U[[1 ;; Length[U] ;; 2]];
(* value of first derivative *)
P = U[[2 ;; Length[U] ;; 2]];
(* initialitazion *)
EQ = ConstantArray[0, 2 n];
(* starting loop *)
Do[
jr = 2 j - 1;
(* values for the current element *)
x1 = X[[j]];
x2 = X[[j + 1]];
If[Length[F] == 0,
f1 = F;
f2 = F,
(* if non-constant *)
f1 = F[[j]];
f2 = F[[j + 1]]
];
If[Length[P] == 0,
p1 = P;
p2 = P,
(* if non-constant *)
p1 = P[[j]];
p2 = P[[j + 1]]
];
(* j-th equation *)
EQ[[jr ;; jr + 3]] = EQ[[jr ;; jr + 3]] +
Flatten[NIntegrate[fun[x, x1, x2, f1, p1, f2, p2], {x, x1, x2}]],
{j, 1, n - 1}
];
(* assembly of system *)
R = EQ - RHS;
(* boundary conditions *)
R[] = U[] - wf0;
R[] = U[] - wp0;
R[[2 n - 1]] = U[[2 n - 1]] - wf1;
R[[2 n]] = U[[2 n]] - wp1;
(*output*)
R
];

where I assumed that $w(x)$ is the exact solution, thus wf0 = w(x0), wp0 = w'(x0), wf1 = w(x1)$and wp1 = w'(x1). Then: (* forcing term *) f1[x_] = w''''[x]; (* integration of forcing term *) F1 = funFF[XX, f1]; (* function to be integrated*) fun1[x_, x1_, x2_, f1_, p1_, f2_, p2_] := Ne2[x, x1, x2, f1, p1, f2, p2] NN2[x, x1, x2]; (* solution *) U1 = UU /. Flatten[FindRoot[SYS1[XX, UU, F1, fun1] == 0, {UU, U0}] where: • XX is the grid • U0 = {\alpha_1, \beta_1, \alpha_2, \beta_2, ... } is the guess • Ne2 contains the interpolation inside an element • NN2 contains the four shape functions • fun is the integrand for each element, i.e$u'' \psi ''\$

This simple example works, but it is really slow even if it's a linear case. How can I improve/optimize the code? I'm new in this kind of programming, so I'm sure I'm missing some speeding-up tricks.

Any suggestion is very appreciated.

Best,

Petrus

• In the finite element method the numerical integration over the element is typically done using very minimal sampling (i.e. one or two quadrature points). You should not be using NIntegrate , rather work out by hand the 1- or 2- point gauss integration formula and use that. – george2079 Jul 11 '14 at 20:39
• Having taken @george2079's suggestion, your code should be ideal for Compile, which will speed it up significantly. – Oleksandr R. Dec 26 '14 at 14:16
• @Petrus, where does your project stand? Do you now have a reliable, reasonably fast procedure? – bbgodfrey Mar 1 '15 at 15:12