Burgers equation

$$ \begin{array} UU_t + UU_x - \nu U_{xx} = 0 \\ U(x,t) \mid_{x=-\infty} = U_1 \\ U(x,t) \mid_{x=+\infty} = U_2 \\ U(x,t) \Bigr|_{t=0} = \left\{ \begin{array}{rcl} U_1 & \mbox{for} & x < a \\ (x - b) \frac{U_1 - U_2}{b-a} +U_2 & \mbox{for} & a < x < b \\ U_2 & \mbox{for} & x > b \end{array}\right. \end{array}$$

Solving with stiff method of lines and Rosenbrock scheme with complex coefficient

Quasi-uniformal mesh will be introduced later

$$ \dfrac{\partial U}{\partial t} = \nu \dfrac{\partial^2 U}{\partial t^2} - U \dfrac{\partial U}{\partial x} $$ $$ \dfrac{dU}{dt} = \nu \left[ \dfrac{U_{i+1} - U_{i}}{x_{i+1} - x_{i}} - \dfrac{U_{i} - U_{i-1}}{x_{i} - x_{i-1}} \right] \frac{2}{x_{i+1} - x_{i-1}} - U_i \dfrac{U_{i+1} - U_{i-1}}{x_{i+1} - x_{i-1}} $$


In [1]:
using Plots
using DifferentialEquations
pyplot()


Out[1]:
Plots.PyPlotBackend()

In [2]:
(a,b,Uₗ,Uᵣ,ν) = (-2.0,2.0,4.0,2.0,3.0)
k = (Uᵣ - Uₗ)/(b-a)
function U₀(x::Float64)
    (x < a) ? Uₗ :
    (x > b) ? Uᵣ :
    k*(x - b) + Uᵣ
end
# Initial condition for U
# plot(map(x -> U₀(x),-5:0.5:5),title="U(x,t=0)",ylims=(-5,5),leg=false)


Out[2]:
U₀ (generic function with 1 method)

constructing quasi-uniformal mesh

$$ \xi = \alpha \tan{\frac{\pi t}{2}} $$

In [3]:
# правильно ли я ввел сетку ?

N,α = (200,3.0)
quasi_uniformal_step(i::Int64) = α*tan(π*(i)/(N+N/4)) #i*6/N #
# x_mesh should be N+2 size for proper end of segments handeling
# x_mesh[1] and [N+1] is ±∞
#                   -500      500
x_mesh = [quasi_uniformal_step(-div(N,2)+i) for i=-1:N]
println("x_mesh length " ,length(x_mesh))
plot(x_mesh[1:N+2],leg=false,title="quasi-uniformal tan() mesh",xlabel="i",ylabel="x[i]") # exept the biggest values at the end


x_mesh length 202
Out[3]:
sys:1: MatplotlibDeprecationWarning: The set_axis_bgcolor function was deprecated in version 2.0. Use set_facecolor instead.

In [4]:
# δx is left difference between adjecent mesh points. U[i+1]-U[i] should correspond to δx[i+1]
δx = [x_mesh[i+1] - x_mesh[i] for i = 1:N]
print("δx size: ",length(δx), " values at the ends: ",δx[1],δx[2]," ",δx[N-1],δx[N])
plot(δx,ylims=(extrema(δx)))


δx size: 200 values at the ends: 0.410695624376316460.3801093839655749 0.32865595167551830.35292476448493204
Out[4]:

In [5]:
# U_inital should be generated not at ±∞ (2:N) while mesh (1:N+1)
U_init = [U₀(x_mesh[i+1]) for i=1:N]
plot(U_init,title="U_init at quasi-uniformal mesh",xlabel="i",ylabel="U₀",ylims=(0,Uₗ+1))


Out[5]:

In [6]:
# quasi-uniformal mesh right part
function rp(t,u,du)
    du[1] = ν * (2/(δx[2]+δx[1]))*( ((u[2] - u[1])/(δx[2])) - ((u[1] - Uₗ)/(δx[1])) )
    - u[1] * (u[2] - Uₗ)/(δx[2]+δx[1])
    for i = 2:N-1 
        du[i] = ν * (2/(δx[i+1]+δx[i]))*( ((u[i+1] - u[i])/(δx[i+1])) - ((u[i] - u[i-1])/(δx[i])) ) 
        - u[i] * (u[i+1] - u[i-1])/(δx[i+1]+δx[i])
    end
    du[N] = ν * (2/(δx[N]+δx[N-1]))*( ((Uᵣ - u[N-1])/(δx[N])) - ((u[N] - u[N-1])/(δx[N-1])) ) 
    - u[N] * (Uᵣ - u[N-1])/(δx[N]+δx[N-1])
end


Out[6]:
rp (generic function with 1 method)

In [7]:
sol = solve(ODEProblem(rp,U_init,(0.0,10.0)));

In [8]:
ts = length(sol[:,1])
println("time steps made: ",ts)
u_result = Array{AbstractVector}(length(sol[:,1]),N)
for i=1:length(sol[:,1])
    u_result[i] = map(tuple,x_mesh[2:N+1],sol[i,:])
end
plot(u_result[ts])


time steps made: 23856
Out[8]:

In [9]:
#anim = Animation();
#for i = 1:length(sol[:,1])
#    p = plot(u_result[i])
#              frame(anim)
#       end
#gif(anim,"qa-un.gif")

In [10]:
h = (x_mesh[N+2]-x_mesh[1])/N
# uniformal mesh right part
function rp2(t,u,du)
    du[1] = ν * (u[2] - 2*u[1] + Uₗ)/h/h - u[1] * (u[2] - Uₗ)/(2*h)
    for i = 2:N-1 
        du[i] = ν * (u[i+1] - 2*u[i] + u[i-1])/h/h - u[i] * (u[i+1] - u[i-1])/(2*h)
    end
    du[N] = ν * (Uᵣ - 2*u[N] + u[N-1])/h/h - u[N] * (Uᵣ - u[N-1])/(2*h)
end


Out[10]:
rp2 (generic function with 1 method)

In [11]:
sol_un = solve(ODEProblem(rp2,U_init,(0.0,10.0)));

In [12]:
ts = length(sol_un[:,1])
println("time steps made: ",ts)
u_result_un = Array{AbstractVector}(ts,N)
for i=1:ts
    u_result_un[i] = map(tuple,x_mesh[2:N+1],sol_un[i,:])
end
plot(u_result_un[ts])


time steps made: 3839
Out[12]:

In [13]:
#anim = Animation();
#for i = 1:length(sol[:,1])
#    p = plot(u_result_un[i])
#              frame(anim)
#       end
#gif(anim,"un.gif")

In [29]:



WARNING: Method definition right_part(Any) in module Main at In[28]:3 overwritten at In[29]:3.
WARNING: Method definition rigth_part_jacobian(Any) in module Main at In[28]:12 overwritten at In[29]:12.
Out[29]:
rigth_part_jacobian (generic function with 1 method)

In [ ]: