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 [ ]: