From the inviscid Burgers' equation: $\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial x} = 0$, I get the discretization $\frac{u_i^{n+1}-u_i^n}{\Delta t}+\frac{(u_i^{n+1})^2-(u_{i-1}^n)^2}{\Delta x}=0$. Having the starting conditions $u_0(x)=sin(\pi x)$, $x\in[0,1]$, I wish to plot $u$ after $n$ timesteps. However, I cannot seem to get it to work and would appreciate any help. Python code:

import numpy as np
import matplotlib.pyplot as plt
import scipy as py
from scipy.sparse.linalg import gmres

x0 = 0 #x min
xe = 1 #x max
N = 100 #Nodes
dx = (xe-x0)/N 
x = py.linspace(x0,xe,N)

t0 = 0 #time
te = 1 #max time
dt = 0.0001 #timestep
M = int((te-t0)/dt) #timesteps

u0 = py.sin(py.pi*x)

ts = 50

def EulerMethod(ustart):
    U = np.zeros((N,M))
    U[:,0] = ustart
    for i in range(N-1): #one step Euler forward
        U[0,1] = 0
        U[i+1,1] = U[i+1,0] - (dt/dx)*(U[i+1,0]**2-U[i,0]**2)
    for n in range(ts):
         U[:,n+2] = NewtonMethod(U[:,n+1],U[:,n]) #Eulers implicit method with Newtons Method        
    return U

def NewtonMethod(u,c):
    #c = u_n
    for i in range(3):
        F = np.zeros((N,1))
        J = np.zeros((N,N))
        J[N-1,N-1] = 1/dt + u[N-1]/dx
        F[0] = (u[0] - c[0])/dt + (u[0]**2)/2*dx
        for j in range(N-1):
            J[j,j] = 1/dt + u[j]/dx
            J[j+1,j] = 1/dt + u[j]/dx
            F[j+1] = (u[j+1] - c[j+1])/dt + (u[j+1]**2 - u[j]**2)/2*dx
        du = gmres(J,-F)[0]
        u = u + du
    return u

UM = EulerMethod(u0)
  • $\begingroup$ I would suggest writing out the math of what you are trying to accomplish so that we can read it. Because from a glance at your code it looks a lot like you're trying to run both forward and implicit Euler at the same time, but your forward Euler is also implemented incorrectly as far as I can tell. $\endgroup$
    – EMP
    Apr 14 '19 at 14:02
  • $\begingroup$ I will try to reformulate it better. Forward Euler is just supposed to be done the first step and I have fixed that now. Thanks for your feedback. $\endgroup$
    – Mathias
    Apr 14 '19 at 16:54
  • $\begingroup$ does it work now? $\endgroup$
    – EMP
    Apr 14 '19 at 18:30
  • $\begingroup$ As something of a side note, your space discretisation in your first sentence isn't a consistent one. $(a^2_i -a^2_{i-1})/\Delta x$ is a one sided approximation of $\frac{d a^2}{dx} = 2a\frac{da}{dx} $. I also wouldn't call the time discretisation you've written implicit Euler, but rather some sort of time splitting method. $\endgroup$
    – origimbo
    Apr 14 '19 at 18:50
  • $\begingroup$ Yea, you also do not need to initialize with a forward Euler step if you're doing 1st order backward Euler. $\endgroup$
    – EMP
    Apr 15 '19 at 0:07

Your Answer

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

Browse other questions tagged or ask your own question.