Implementation of Backward-Euler scheme, Newton-Raphson iteration scheme to time dependent nonlinear differential equation

Question

I just knew how to do Newton-Raphson iteration in time-independent 1D nonlinear differential equation. Then I applied to time-dependent 1D nonlinear differential equation, and I got confused.

Below is just test time-dependent 1D nonlinear differential equation I want to solve (I just made this up from heat equation, to make it nonlinear) $$ \frac{\partial u}{\partial t}+\frac{\partial^2 u}{\partial x^2}+u^2=0 $$ For this equation, I assumed I would need initial condition, and boundary condition $$ u(x,t=0)=\left\{ \begin{aligned} 1\qquad&(-1\leq x\leq 1)\\ 0\qquad&\text{otherwise} \end{aligned}\right.,\qquad u(-5,t)=0,\quad u(5,t)=0 $$ I set the $x$ range of $-5 \leq x \leq 5$, and planned to create code in Matlab.

I started with space discretization by using second-order difference $$ \frac{\partial u}{\partial t}+\frac{u_{i+1}-2u_{i}+u_{i-1}}{\Delta x^2}+u_{i}^2=0 $$ where $i$ denotes space discretization number

Then I applied Backward-Euler discretization for time discretization, which is $$ \frac{d q(x,t)}{dt}=f(t,x,q(x,t))\qquad\longrightarrow\qquad \frac{q_{j+1}-q_{j}}{\Delta t}=f_{j+1} $$ where $j$ denotes time discretization number (this is just intermediate step for me. I am planning to apply Trapezoidal Rule - Second order Backward Difference Formula (TR-BDF2) later) Now applying this, the equation looks like $$ \frac{u_{i,j+1}-u_{i,j}}{\Delta t}+\frac{u_{i+1,j+1}-2u_{i,j+1}+u_{i-1,j+1}}{\Delta x^2}+u_{i,j+1}^2=0 $$ Again, $i$ and $ j$ denotes space and time discretization respectively.

Then I applied Newton-Raphson(NR) scheme, which is for give equation $$ F(u)=0 $$ the solution is determined by iteration $$ F'(u^{k})\delta u^{k}=-F(u^{k}), \qquad \qquad u^{k+1}=u^{k}+\delta u^{k} $$ (I'm running out of space for discretization) where k denotes NR iteration number, and if $F(u)$ was system of equations, $F'(u)\delta u$ should be Jacobian. Do I have to set initial and boundary condition for $\delta u_{i,j}$ as well?

Now my equation looks like $$ \frac{\delta u_{i,j+1}^{k}-\delta u_{i,j}^{k}}{\Delta t} + \frac{\delta u_{i+1,j+1}^{k}-2\delta u_{i,j+1}^{k}+\delta u_{i-1,j+1}^{k}}{\Delta x^2} +2u_{i,j+1}^{k}\delta u_{i,j+1}^{k} = -\frac{u_{i,j+1}^{k}-u_{i,j}^{k}}{\Delta t} -\frac{u_{i+1,j+1}^{k}-2u_{i,j+1}^{k}+u_{i-1,j+1}^{k}}{\Delta x^2} -(u_{i,j+1}^{k})^2 $$

This is where I got stuck.

How do I proceed after this?

Am I missing something? like another boundary/initial conditions?

Because I only know $u_{i,1}^{0}$ and $\delta u_{i,1}^{0}$.

There are like, 11 unknowns due to backward differences ($j+1$).

Some literature says that I should solve this equation for "each time step"

Or, is this right way to apply Backward-Euler scheme and NR scheme to the time-dependent nonlinear differential equation?

Is there any good example solving time dependent nonlinear differential equation with Newton-Raphson iteration?

Answer 1

You seem to be confused with what are unknowns in your algebraic equations (as mentioned also in the comment by origimbo).

Let say you have $N$ inner nodes for space discretization, i.e. $\Delta x = 10/(N+1)$. The values $u_{i,1}$ shall be obtained from initial condition (a more usual notation is $u_{i,0}$).

The notion "you have to solve your algebraic equations in each time step" means that in each time step the values $u_{i,j}$ are known (either from initial condition or from previous time step) and your unknowns (the values to be found) are only $u_{i,j+1}$. Now this is the task for NR method.

When applying NR method, you have to choose the values $u^0_{i,j+1}$, typically you take $u^0_{i,j+1}=u_{i,j}$. It means $F(u)=0$ represents $N$ nonlinear algebraic equations having $N$ unknowns $u_{1,j+1}$ up to $u_{N,j+1}$ having the form $$ \frac{\delta u_{i,j+1}^{k}}{\Delta t} + \frac{\delta u_{i+1,j+1}^{k}-2\delta u_{i,j+1}^{k}+\delta u_{i-1,j+1}^{k}}{\Delta x^2} +2u_{i,j+1}^{k}\delta u_{i,j+1}^{k} = -\frac{u_{i,j+1}^{k}-u_{i,j}}{\Delta t} -\frac{u_{i+1,j+1}^{k}-2u_{i,j+1}^{k}+u_{i-1,j+1}^{k}}{\Delta x^2} -(u_{i,j+1}^{k})^2 $$

In above, the values $u_{i,j}$, $u_{*,j+1}^k$ are known, the unkowns in above are only three (but in each equation different ones), namely $\delta u_{i-1,j+1}^k$, $\delta u_{i,j+1}^k$, and $\delta u_{i+1,j+1}^k$. The boundary conditions are $\delta u_{0,j+1}^k=0$ and $\delta u_{N+1,j+1}^k=0$.

By the way, the sign in your PDE before the diffusion term is wrong, it should be minus.