Quasi-uniformal mesh will be introduced later SMOL: $$ \dfrac{\partial U}{\partial t} = \nu \dfrac{\partial^2 U}{\partial x^2} - U \dfrac{\partial U}{\partial x} $$ $$ \left\{ \begin{array}{lc1} \dfrac{dU_i}{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}} \quad i=\overline{1,N} \\ U_0(t) = U_l \\ U_{N+1}(t) = U_r \\ U_i(0) = U_0(x_i)\end{array}\right.$$

$$ \dfrac{dU_i}{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}} = F_i$$$$ \xi_i = tan(\frac{i \pi}{N}) \quad i=-\frac{N}{2},\frac{N}{2} $$$$ \delta x_i = \xi_{i+1} - \xi_i \quad \mbox{for} \quad i= \overline{1,N} $$

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

```
Out[1]:
```

```
In [2]:
```(a,b,Uₗ,Uᵣ,ν) = (0,π,2.,0.,30.) #inital params
N,α = (200,3.0) #quasi-uniformal mesh
k = (Uᵣ - Uₗ)/(b-a)
tanh_i = false
function U₀(x::Float64)
(x < a) ? Uₗ :
(x > b) ? Uᵣ :
1 + cos(x)
end
#ξ(i) = 1*(-N/2 + i)/N + 0.5
#if tanh_i
# function U₀(x::Float64)
# - (((x)^2 - x -2 ) - 6 * tanh( -3 * ( (x-0.6))))
# end
# ξ(i) = 1*(-N/2 + i)/N + 0.5
# space_mesh = [ξ(i) for i=0:N+1]
# Uₗ = U₀(space_mesh[1])
# Uᵣ = U₀(space_mesh[N+2])
#end -N/2.0 +
ξ(i) = 3*π*(-N/2 + i)/N
#ξ(i) = α*tan(π*(-N/2 + i)/(N + N/64))
space_mesh = [ξ(i) for i=0:N+1];

```
In [3]:
```δx = [space_mesh[i+1] - space_mesh[i] for i = 1:N]
U_init = [U₀(space_mesh[i+1]) for i=1:N]
pui = plot(U_init,title="U_init at quasi-uniformal mesh",xlabel="i",ylabel="U₀",ylims=(0,max(Uᵣ,Uₗ)+1));
pξ = plot(space_mesh,ylabel="ξ(i)",xlabel="i",leg=false);
pδ = plot(δx,ylabel="δx(i)",xlabel="i",leg=false,ylims=(extrema(δx)));
println("Uᵢ ",length(U_init)," space_mesh ",length(space_mesh)," δx ",length(δx))
println("ξ ",space_mesh)
println("δx",δx)
plot(pξ,pδ,pui)

```
Out[3]:
```

$$ \dfrac{dU_i}{dt} = \nu \left[ \dfrac{U_{i+1} - U_{i}}{\delta x_{i+1}} - \dfrac{U_{i} - U_{i-1}}{\delta x_{i}} \right] \frac{2}{\delta x_{i+1} + \delta x_{i}} - U_i \dfrac{U_{i+1} - U_{i-1}}{\delta x_{i+1}+\delta x_{i}} = F_i$$$$ \dfrac{\partial F_1}{\partial x_1} = \frac{2\nu}{\delta x_{2}+\delta x_1} \left( -\frac{1}{\delta x_2} -\frac{1}{\delta x_1} \right) - \frac{U_2 - U_l}{\delta x_2 + \delta x_1}$$$$ \dfrac{\partial F_1}{\partial x_2} = \frac{2\nu}{\delta x_{2}+\delta x_1} \left( \frac{1}{\delta x_2} \right) - \frac{U_1}{\delta x_2 + \delta x_1}$$$$ \dfrac{\partial F_i}{\partial x_{i-1}} = \frac{2\nu}{\delta x_{i+1}+\delta x_i} \left( \frac{1}{\delta x_i} \right) + \frac{U_i}{\delta x_{i+1}+\delta x_i} $$$$ \dfrac{\partial F_i}{\partial x_i} = \frac{2\nu}{\delta x_{i+1}+\delta x_i} \left( -\frac{1}{\delta x_{i+1}} -\frac{1}{\delta x_{i}} \right) - \frac{U_{i+1} - U_{i-1}}{\delta x_{i+1}+\delta x_i}$$$$ \dfrac{\partial F_i}{\partial x_{i+1}} = \frac{2\nu}{\delta x_{i+1}+\delta x_i} \left( \frac{1}{\delta x_{i+1}} \right) - \frac{U_i}{\delta x_{i+1}+\delta x_i}$$$$ \dfrac{\partial F_N}{\partial x_{N-1}} = \frac{2\nu}{\delta x_{i+1}+\delta x_i} \left( \frac{1}{\delta x_{N-1}} \right) + \frac{U_N}{\delta x_N + \delta x_{N-1}} $$$$ \dfrac{\partial F_N}{\partial x_{N}} = \frac{2\nu}{\delta x_{i+1}+\delta x_i} \left( -\frac{1}{\delta x_N} -\frac{1}{\delta x_{N-1}} \right) - \frac{U_r - U_N}{\delta x_N + \delta x_{N-1}}$$

```
In [4]:
```function right_part(U::Vector)
r = similar(U)
r[1] = ν*(2.0/(δ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
r[i] = ν*(2.0/(δ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
r[N] = ν*(2.0/(δx[N]+δx[N-1]))*((Uᵣ-U[N])/(δx[N]) - (U[N] - U[N-1])/(δx[N-1])) - U[N] * (Uᵣ-U[N-1])/(δx[N]+δx[N-1])
return r
end
function right_part_jacobian(U)
dl = Array{Float64}(N-1)
du = Array{Float64}(N-1)
d = Array{Float64}(N)
du[1] = ν*(2/(δx[2]+δx[1]))*(-1/δx[2]-1/δx[1]) - (U[2] - Uₗ)/(δx[2]+δx[1])
d[1] = ν*(2/(δx[2]+δx[1]))*(1/δx[2]) - U[1]/(δx[2]+δx[1])
for i = 2:N-1
du[i] = ν*(2/(δx[i+1]+δx[i]))*(1/δx[i+1]) + U[i]/(δx[i+1]+δx[i])
d[i] = ν*(2/(δx[i+1]+δx[i]))*(-1/δx[i+1]-1/δx[i]) - (U[i+1] - U[i-1])/(δx[i+1]+δx[i])
dl[i] = ν*(2/(δx[i+1]+δx[i]))*(1/δx[i+1]) - U[i]/(δx[i+1]+δx[i])
end
dl[N-1] = ν*(2/(δx[N]+δx[N-1]))*(1/δx[N-1]) + U[N]/(δx[N]+δx[N-1])
d[N] = ν*(2/(δx[N]+δx[N-1]))*(-1/δx[N]-1/δx[N-1]) - (Uᵣ - U[N])/(δx[N]+δx[N-1])
return Tridiagonal(dl,d,du)
end
g = x -> ForwardDiff.jacobian(right_part, x)

```
Out[4]:
```

```
In [5]:
```solution = U_init
#solution[N] += (Uᵣ - solution[N-1])/2
println("Uₗ Uᵣ ",Uₗ ," ", Uᵣ)
println("U_initial: ",solution)
println("F_i = rigth_part(U_initial): ",right_part(solution))
println(real((eye(N) - complex(eye(N),eye(N)).*(0.001/2))*right_part_jacobian(solution))\ right_part(solution))
# println(Uᵣ," ",solution[N]," ",solution[N-1])
# Некорректная из-за близости значений в x=a, х=b и соседних точках
plot(right_part(solution),lab="F(U_initial)",title="incorrect derivative approximation \n of partialy continious function")

```
Out[5]:
```

```
In [6]:
```tspan = (0.,1.,1000.,0.); t=0
τ = (tspan[2] - tspan[1])/tspan[3]
solution = U_init
anim = Animation()
while t<=tspan[3]
if (t%50==0)
plot(solution,ylims=(-10.,10.))
frame(anim)
end
W = ((eye(N) - complex(eye(N),eye(N)).*(τ/2))*right_part_jacobian(solution))\ right_part(solution)
solution += τ * real(W)
solution[1] = solution[2] = solution[3] = Uₗ
solution[end] = solution[end-1] = solution[end-2] = Uᵣ
t+=1
end

```
In [7]:
```plot(solution)

```
Out[7]:
```

```
In [8]:
```gif(anim,"manual_burgers.gif")

```
Out[8]:
```

```
In [9]:
``````
solution
```

```
Out[9]:
```

```
In [10]:
```function rp(t,u,du)
r = Vector(N)
du[1] = ν*(2.0/(δ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.0/(δ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.0/(δx[N]+δx[N-1]))*((Uᵣ-u[N])/(δx[N]) - (u[N] - u[N-1])/(δx[N-1])) - u[N] * (Uᵣ-u[N-1])/(δx[N]+δx[N-1])
end
sol = solve(ODEProblem(rp,U_init,(0.,1.)),ImplicitEuler(),dt=0.01)

```
Out[10]:
```

```
In [11]:
```println(indices(sol[:,1]))
plot(sol[end,:])

```
Out[11]:
```

```
In [12]:
```sol[end,:]

```
Out[12]:
```

```
In [ ]:
```