RK4 gives nan for finite element galerkin method with 8+ basis

Question

Exact Solution = $e^t \sin(\pi(x)) $

$f(x,t) = e^t(1-\beta \pi^2)\sin(\pi(x))$

$H = 1/N$

$ \phi_0 = (H - x)/H \hspace{5mm} in \hspace{5mm} [0, H] \hspace{5mm} else \hspace{5mm} 0 $

$ \phi_N = (x - (N-1)H )/H \hspace{5mm} in \hspace{5mm} [(N-1)H, NH] \hspace{5mm} else \hspace{5mm} 0 $

$ \phi_k = (x -(k - 1)H )/H \hspace{5mm} in \hspace{5mm} [kH, (k+1)H],\hspace{5mm} ((k + 1)H - x )/H \hspace{5mm} in \hspace{5mm} [kH, (k+1)H] \hspace{5mm} else \hspace{5mm} 0 $

I am trying to solve the equation $ u_t + \beta u_{xx} = f(x,t)$. I have used $ u = \sum_{0}^{N} a_i(t) \phi_i(x) $ as the approximate solution where $\phi_i$ is defined by the triangular hat functions in [0, 1]. I 'm using finite element analysis and using the sense of weak derivatives I have calculated the various stiffness matrices that would be required. I have used hat functions as the test functions too.

After a bit of solving I let 2 tridiagonal matrices and one other with 4 entries. Now I have to solve the system $$ \alpha^{'}_{(n \times 1)}= A^{-1}_{(n \times n)}({D_{(n \times n)}.y_{(n \times 1)} + e^t F_{(n \times 1)}})$$ where $'$ denotes derivatives and subscripts the dimensions of the matrices and $A$, $D$, $F$ are matrices.

I had been trying to solve this system by

• (1) RK4 (self implemented),
• (2) ODE45 (scipy's RK45 and LSODA),
• (3) Backward Euler (self implemented).

For (1) the solution blows up for n=8 basis or elements and further. For (2) LSODA I get the error initially decreasing but increasing after 100 or so basis and overall the $\log (error)$ vs. $\log(H)$ plot is increasing(unexpected). For RK45 I get nan values both on my solver and scipy's too. For (3) I get decreasing error plot for $\log(error)$ vs $\log(H)$.

Main problem is why the Runge-Kutta method. It is supposed to be the best. I was taking time step as $H^2$ in Runge-Kutta. Distance step is H as N + 1 points line on the x -axis with height of triangle 1 .

Please help me out. I can share the code and whatever's required to elaborate further.

Answer 1

There were some clarifications reached in the comments, esp. establishing that $β=-1$ is negative and including the actual code, that give a different weight to the points raised in the previous answer.

Construction of the test case

This experiment in the accuracy of the method of lines is based on the heat equation operator $L[u]=u_t-u_{xx}$ with boundary condition operator $R[u]=(u(0),u(1))$.

Method of Manufactured Solutions (MMS)

In this method one constructs a test problem as $L[u]=f=L[p]$, $R[u]=R[p]$, $u|_{t=0}=p|_{t=0}$, where $p$ is some smooth function that serves as the known exact solution. Here

1. $p(x,t)=e^t\sin(\pi x)$ giving homogeneous boundary conditions and
2. $p(x,t)=e^t\sin(x)$ giving an inhomogeneous condition at the $x=1$ boundary

By using second order accurate approximation of the space derivatives, the expected error is $O(\Delta t^p+\Delta x^2)$ where $p$ is the order of some fixed-step ODE solver.

The space discretization error can be removed by applying the MMS approach to the discretization $L_h$ of the differential operator, with $\Delta x=h$. Then in $L_h[u_h]=f_h=L_h[p_h]$, $R[u_h]=R[p_h]$ etc. the error of $u_h$ against $p_h$ is only depending on the time. $p_h$ is some preliminary space discretization of $p$, for instance it could be the projection on the piecewise linear functions.

Finite element method

The weak solution of the problem satisfies $$ \partial_t\int_0^1 uw\,dx+\int_0^1u_xw_x\,dx=\int_0^1fw $$ for all differentiable functions $w$ with $w(x=0)=w(x=1)=0$. The finite element solution has $u$ and $w$ in the space of piecewise linear functions in $x$ direction, over the uniform grid of step size $\Delta x=h=\frac1n$, $x_k=kh$.

Set $u(x,t)=\sum_{k=0}^N u_k(t)\phi_k(x)$, $\phi_k(x)=\max(0, 1-|x/h-k|)$ the hat functions. This results in the well-known coefficient sequences $\frac16[1,4,1]$ and $\frac1{h^2}[-1,2,-1]$ for the scalar products of the $\phi_k$ and the $\partial_x\phi_k$. For $k=1,...,N-1$ this thus results in equations $$ \frac{u_{k+1}'(t)+4u_k'(t)+u_{k-1}'(t)}6 + \frac{-u_{k+1}(t)+2u_k(t)-u_{k-1}(t)}{h^2}=\frac1h\int_0^1f(x,t)\phi_k(x)\,dx=F_k(t). \tag{FEM} $$ The integral on the right evaluates for $p(x,t)=e^t\sin(wx)\implies f(x,t)=e^t(1+w^2)\sin(wx)$ to $$ F_k(t)=e^t(1+w^2)\left(\frac{\sin(wh/2)}{wh/2}\right)^2\sin(wx_k). $$ Also note that $a\sin(X+H)+b\sin(X)+a\sin(X-H)=(b+2a\cos(H))\sin(X)$.

Homogeneous case

The continuous problem, the exact solution, the discretization and all difference and differential operators are symmetric so that the solution will always be a multiple of $\sin(\pi x)$ up to floating point errors that are continuously dampened by the properties of the heat equation. The exact solution of the discretized equation will thus have the form $u_k(t)=a(t)\sin(\pi x_k)$, $a(0)=1$, and satisfy the differential equation $$ \frac{(4+2\cos(\pi h))}6a'(t)\sin(\pi x_k) +\frac{2(1-\cos(\pi h))}{h^2}a(t)\sin(\pi x_k) =e^t(1+\pi^2)\left(\frac{\sin(\pi h/2)}{\pi h/2}\right)^2\sin(\pi x_k) \\\iff\\ a'(t) + \pi^2 c(h)a(t) = c(h)(1+\pi^2) e^t,~~ c(h)=\frac{\left(\frac{\sin(\pi h/2)}{\pi h/2}\right)^2}{\frac{(4+2\cos(\pi h))}6}=1+O(h^2) $$ The numerical solver for the system will effectively solve this scalar equation with solution $$ a(t)=\frac{1+\pi^2}{1+c(h)\pi^2} (e^t-e^{-c(h)\pi^2 t})+e^{-c(h)\pi^2 t} $$ The error against the prescribed function $e^t$ is proportional to $c(h)-1=O(h^2)$ which is also observed in the numerical experiments.

Inhomogeneous case

In the equation (FEM) for the index $k=n-1$ the outer node $u_n(t)$ is not zero and also not constant, so it has a contribution to both the time derivative term and to the second order space difference term. As $u_n(t)$ is not a variable in the state vector for the method-of-lines ODE system, these components have to be added to the forcing term on the right side. In code this could look like

x = np.linspace(0, 1, N+1); #0 - 1
H = x[1]-x[0];
S = np.sin(x[1:-1]);
p = lambda t:math.exp(t) * self.alpha
b_c_start = lambda t : 0
b_c_end = lambda t :math.exp(t)*math.sin(1)
                    
A = sp.diags([1/6, 2/3, 1/6], [-1, 0, 1], shape=(N-1, N-1));
Ainv = sp.linalg.splu;
D = sp.diags([-1, 2, -1], [-1, 0, 1], shape=(N-1, N-1));
eN = np.zeros(N-1); eN[N]=1; 
DN=-eN;
AN = eN/6;
F = (math.sin(H/2)/(H/2))**2 * (1+1) * S #incomplete without exp(t)

def MoLfunc (self, t,  y):
    # vectorized = False, default
    yN = b_c_end(t)
    Dy = H**-2*(D.dot(y) + DN*yN);
    return Ainv.solve(np.exp(t)*F - AN*yN - Dy);  

Also note that the error has to be computed as of the continuous function over $[0,1]$. This means that the Euclidean norm of the difference of the solution vector to the prescriped function values has to be scaled by $\sqrt{h}$ to approximate the $L^2$ norm.

Conclusion

With all these corrections the integration proceeds without problem and the estimated exponent of $h$ in the error is close to 2, $1.996$ in the homogeneous case and $1.99986$ in the inhomogeneous case. The differences are due to details in the computations , the ODE solver used and the error tolerances passed to it.