# Solving an iterative, implicit Euler method in MATLAB

I'm trying to solve an iterative problem that includes an implicit (backwards) Euler method to find successive time values for a given function. The numerical problem is shown here:

\begin{aligned} \frac{u_{r,s+1}-u_{r,s}}{\Delta t}&+u_{r,s}\frac{u_{r+1,s+1}-u_{r-1,s+1}+u_{r+1,s}-u_{r-1,s}}{4\Delta x}\\ &\quad+\frac{\eta^*_{r+1,s+1}-\eta^*_{r-1,s+1}+\eta_{r+1,s}-\eta_{r-1,s}}{4\Delta x}\\ &=\frac{1}{3}(\alpha x)^2\frac{u_{r+1,s+1}-2u_{r,s+1}+u_{r-1,s+1}-u_{r+1,s}+2u_{r,s}-u_{r-1,s}}{\Delta x^2 \Delta t}\\ &\quad+\alpha^2x\frac{u_{r+1,s+1}-u_{r-1,s+1}-u_{r+1,s}+u_{r-1,s}}{\Delta x\Delta t}, \end{aligned}

where $s$ is the time iteration and $r$ is the location iteration, $x$ is the true location ($x=r\,\Delta x$), and $\alpha$, $\Delta x$, and $\Delta t$ are constants. The $\eta$ values are already solved, so I need to solve for $u(r,s+1)$. However, the solution includes values for $u(r-1,s+1)$ and $u(r+1,s+1)$, which requires this to be solved as a system of equations over $r$ for each $s$. I have boundary conditions of $u$ for all $r$ when $s=1$, as well as u at the first and last $r$ values for all $s$ (i.e. $u(1,s)$ and $u(r_\max,s)$). I'm not sure what to do from here though to solve for the "interior" $u$ values at $s+1$.

For reference, the original, non-discretized equation is here:

$$\frac{\partial\overline{u}}{\partial t}+\overline{u}\frac{\partial\overline{u}}{\partial x}+\frac{\partial\eta}{\partial x}=\frac{1}{3}(\alpha x)^2 \frac{\partial^3\overline{u}}{\partial x^2\partial t}+\alpha^2 x\frac{\partial^2 \overline{u}}{\partial x \partial t}.$$

• Move all s+1 terms to the left, all others to the right, develop a linear system A*u=b for the equations, change the first and last equations (rows) to match your boundary conditions, and solve: u = A\b. – TroyHaskin Apr 14 '17 at 1:34
• Hi Troy, I don't think I can do that because the second term isn't linear (i.e. it multiplies u(s+1) terms with a u(s) term, so I can't separate all the (s+1) terms from everything else). – Jonathan Lala Apr 14 '17 at 18:01
• The u_{r,s} terms are all known; therefore, they are simply coefficients. Thus, the system is linear in the s+1 terms. (In truth, there is more than likely some linearization at work to make this so since a big problem with momentum equations is the u^2 nonlinearity. But I suspect this was precluded in a low-Mach number sense to make the problem less difficult to solve.) – TroyHaskin Apr 15 '17 at 2:16
• Is there a reason you want to use an Euler method? You can probably do a similar discretization and throw it to ode15s. – Chris Rackauckas Apr 15 '17 at 18:53
• Why does your $x$ derivative involve values at two different times? – nicoguaro Apr 16 '17 at 0:58

Thanks to TroyHaskin, I've realized that this is in reality a linear problem, since the u(s) values are known from the initial conditions, and is therefore easy to solve in Matlab. Thanks to everyone else for help as well.

EDIT: To be more specific, the system can be solved linearly by separating the u(s+1) terms and their coefficients from everything else. The solution takes the form [Aw,Ap,Ae]u = Q, where u = [u(r-1,s+1),u(r,s+1),u(r+1,s+1)]^T. Because this is a tridiagonal matrix, it can be solved with minimum storage (note that I have changed location r to i and time s to n in my code, just for personal consistency):

    A_E = zeros(N,1);
A_W = zeros(N,1);
A_P = zeros(N,1);
Q = zeros(N+1,1);
Q(1) = u(1,n+1);
A_P(1) = 1;

for i = 2:N
A_W(i) = (-u(i,n)/(4*dx))-(((alpha*x(i))^2)/(3*(dx^2)*dt))+((alpha^2)*x(i)/(dx*dt));
A_P(i) = (1/dt)+((2*((alpha*x(i))^2))/(3*(dx^2)*dt));
A_E(i) = (u(i,n)/(4*dx)) - (((alpha*x(i))^2)/(3*(dx^2)*dt)) - ((alpha^2)*x(i)/(dx*dt));
Q(i) = u(i,n)/dt - (u(i,n)/(4*dx))*(u(i+1,n)-u(i-1,n)) ...
- (eta(i+1,n+1)-eta(i-1,n+1)+eta(i+1,n)-eta(i-1,n))/(4*dx) ...
+ (((alpha*x(i))^2)/(3*(dx^2)*dt))*(-u(i+1,n)+(2*u(i,n))-u(i-1,n)) ...
+ ((alpha^2)*x(i)/(dx*dt))*(-u(i+1,n)+u(i-1,n));
end

%Gaussian Elimination
for i = 2:N
A_P(i) = A_P(i) - (A_W(i)*A_E(i-1)/A_P(i-1));
Q(i) = Q(i) - (A_W(i)*Q(i-1)/A_P(i-1));
end

%Back Substitution
for i = N:-1:2
u(i,n+1) = (Q(i) - (A_E(i)*u(i+1,n+1)))/A_P(i);
end

• Just a friendly heads up: you should post your solution since this is currently not an answer, and that is frowned upon. – TroyHaskin Apr 19 '17 at 6:39
• Fixed. Thanks for the heads up; I'm new to stack exchange (and compsci in general). – Jonathan Lala Apr 22 '17 at 12:35