2
$\begingroup$

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

$\endgroup$

migrated from stackoverflow.com Apr 15 '17 at 15:28

This question came from our site for professional and enthusiast programmers.

  • 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 '17 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 '17 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 '17 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$ – Chris Rackauckas Apr 15 '17 at 18:53
  • 1
    $\begingroup$ Why does your $x$ derivative involve values at two different times? $\endgroup$ – nicoguaro Apr 16 '17 at 0:58
2
$\begingroup$

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
$\endgroup$
  • 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 '17 at 6:39
  • $\begingroup$ Fixed. Thanks for the heads up; I'm new to stack exchange (and compsci in general). $\endgroup$ – Jonathan Lala Apr 22 '17 at 12:35

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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