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}. $$

  • 3
    $\begingroup$ 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. $\endgroup$
    – TroyHaskin
    Apr 14, 2017 at 1:34
  • $\begingroup$ 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). $\endgroup$
    – Jonathan Lala
    Apr 14, 2017 at 18:01
  • 1
    $\begingroup$ 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.) $\endgroup$
    – TroyHaskin
    Apr 15, 2017 at 2:16
  • 2
    $\begingroup$ Is there a reason you want to use an Euler method? You can probably do a similar discretization and throw it to ode15s. $\endgroup$ Apr 15, 2017 at 18:53
  • 1
    $\begingroup$ Why does your $x$ derivative involve values at two different times? $\endgroup$
    – nicoguaro
    Apr 16, 2017 at 0:58

1 Answer 1


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));

        %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));

        %Back Substitution
        for i = N:-1:2
            u(i,n+1) = (Q(i) - (A_E(i)*u(i+1,n+1)))/A_P(i);
  • 1
    $\begingroup$ Just a friendly heads up: you should post your solution since this is currently not an answer, and that is frowned upon. $\endgroup$
    – TroyHaskin
    Apr 19, 2017 at 6:39
  • $\begingroup$ Fixed. Thanks for the heads up; I'm new to stack exchange (and compsci in general). $\endgroup$ Apr 22, 2017 at 12:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.