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_2 - U_1}{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 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]:





Plots.PyPlotBackend()



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)









    



Uᵢ 200 space_mesh 202 δx 200
ξ [-4.71239,-4.66527,-4.61814,-4.57102,-4.52389,-4.47677,-4.42965,-4.38252,-4.3354,-4.28827,-4.24115,-4.19403,-4.1469,-4.09978,-4.05265,-4.00553,-3.95841,-3.91128,-3.86416,-3.81704,-3.76991,-3.72279,-3.67566,-3.62854,-3.58142,-3.53429,-3.48717,-3.44004,-3.39292,-3.3458,-3.29867,-3.25155,-3.20442,-3.1573,-3.11018,-3.06305,-3.01593,-2.96881,-2.92168,-2.87456,-2.82743,-2.78031,-2.73319,-2.68606,-2.63894,-2.59181,-2.54469,-2.49757,-2.45044,-2.40332,-2.35619,-2.30907,-2.26195,-2.21482,-2.1677,-2.12058,-2.07345,-2.02633,-1.9792,-1.93208,-1.88496,-1.83783,-1.79071,-1.74358,-1.69646,-1.64934,-1.60221,-1.55509,-1.50796,-1.46084,-1.41372,-1.36659,-1.31947,-1.27235,-1.22522,-1.1781,-1.13097,-1.08385,-1.03673,-0.989602,-0.942478,-0.895354,-0.84823,-0.801106,-0.753982,-0.706858,-0.659734,-0.612611,-0.565487,-0.518363,-0.471239,-0.424115,-0.376991,-0.329867,-0.282743,-0.235619,-0.188496,-0.141372,-0.0942478,-0.0471239,0.0,0.0471239,0.0942478,0.141372,0.188496,0.235619,0.282743,0.329867,0.376991,0.424115,0.471239,0.518363,0.565487,0.612611,0.659734,0.706858,0.753982,0.801106,0.84823,0.895354,0.942478,0.989602,1.03673,1.08385,1.13097,1.1781,1.22522,1.27235,1.31947,1.36659,1.41372,1.46084,1.50796,1.55509,1.60221,1.64934,1.69646,1.74358,1.79071,1.83783,1.88496,1.93208,1.9792,2.02633,2.07345,2.12058,2.1677,2.21482,2.26195,2.30907,2.35619,2.40332,2.45044,2.49757,2.54469,2.59181,2.63894,2.68606,2.73319,2.78031,2.82743,2.87456,2.92168,2.96881,3.01593,3.06305,3.11018,3.1573,3.20442,3.25155,3.29867,3.3458,3.39292,3.44004,3.48717,3.53429,3.58142,3.62854,3.67566,3.72279,3.76991,3.81704,3.86416,3.91128,3.95841,4.00553,4.05265,4.09978,4.1469,4.19403,4.24115,4.28827,4.3354,4.38252,4.42965,4.47677,4.52389,4.57102,4.61814,4.66527,4.71239,4.75951]
δx[0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239,0.0471239]






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





(::#5) (generic function with 1 method)



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")









    



Uₗ Uᵣ 2.0 0.0
U_initial: [2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,1.99889,1.99556,1.99002,1.98229,1.97237,1.96029,1.94609,1.92978,1.9114,1.89101,1.86863,1.84433,1.81815,1.79016,1.76041,1.72897,1.69591,1.66131,1.62524,1.58779,1.54902,1.50904,1.46793,1.42578,1.38268,1.33874,1.29404,1.24869,1.20279,1.15643,1.10973,1.06279,1.01571,0.968589,0.921541,0.874667,0.828071,0.781857,0.736127,0.690983,0.646525,0.602852,0.560061,0.518246,0.477501,0.437917,0.39958,0.362576,0.326987,0.292893,0.260369,0.229487,0.200315,0.172919,0.14736,0.123693,0.101972,0.0822454,0.064556,0.0489435,0.0354426,0.0240832,0.0148907,0.0078853,0.00308267,0.00049344,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
F_i = rigth_part(U_initial): [0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,-14.9737,-29.867,-29.6736,-29.4149,-29.0919,-28.7054,-28.2568,-27.7472,-27.178,-26.5508,-25.8671,-25.1286,-24.3373,-23.4949,-22.6035,-21.6651,-20.6819,-19.6561,-18.59,-17.4859,-16.3462,-15.1734,-13.97,-12.7385,-11.4814,-10.2014,-8.90113,-7.5832,-6.2503,-4.90513,-3.55039,-2.1888,-0.823066,0.544074,1.9099,3.2717,4.62675,5.97236,7.30583,8.6245,9.92572,11.2068,12.4653,13.6985,14.9039,16.079,17.2214,18.3287,19.3984,20.4284,21.4163,22.36,23.2573,24.1063,24.905,25.6514,26.3439,26.9807,27.5602,28.0809,28.5415,28.9407,29.2773,29.5504,29.7589,29.9022,28.313,6.66612,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]
[0.0204431,0.0102055,0.020379,0.0305206,0.0406304,0.0507084,0.0607549,0.0707698,0.0807533,0.0907055,0.100627,0.110516,0.120375,0.130203,0.14,0.149767,0.159502,0.169208,0.178882,0.188527,0.198141,0.207725,0.217279,0.226803,0.236297,0.245761,0.255196,0.264601,0.273977,0.283323,0.292639,0.301927,0.311186,0.320415,0.329615,0.338787,0.34793,0.357044,0.36613,0.375187,0.384215,0.393216,0.402188,0.411132,0.420048,0.428936,0.437796,0.446628,0.455433,0.46421,0.472959,0.481681,0.490376,0.499043,0.507683,0.516296,0.524882,0.533441,0.541974,0.550479,0.558958,0.56741,0.575836,0.584235,0.592608,0.600955,0.609276,0.61757,0.625838,0.634081,0.642297,0.650488,0.658653,0.666793,0.674907,0.682995,0.691059,0.699097,0.707109,0.715097,0.723059,0.730997,0.738909,0.746797,0.75466,0.762499,0.770312,0.778102,0.785867,0.793607,0.801323,0.809015,0.816683,0.824327,0.831947,0.839543,0.847115,0.854663,0.862188,0.869689,0.876059,0.880198,0.882126,0.881867,0.879452,0.874913,0.868288,0.859621,0.848959,0.836351,0.821854,0.805525,0.787427,0.767627,0.746193,0.723199,0.698718,0.672831,0.645618,0.617163,0.587552,0.556872,0.525216,0.492673,0.459338,0.425306,0.390673,0.355537,0.319996,0.284149,0.248094,0.211933,0.175765,0.13969,0.103806,0.0682149,0.0330134,-0.00170025,-0.0358292,-0.0692779,-0.101952,-0.133759,-0.164607,-0.194408,-0.223074,-0.25052,-0.276664,-0.301426,-0.324728,-0.346497,-0.36666,-0.38515,-0.401902,-0.416855,-0.429951,-0.441136,-0.45036,-0.457578,-0.462747,-0.465831,-0.466795,-0.465612,-0.462258,-0.456713,-0.448962,-0.438996,-0.426932,-0.414376,-0.401819,-0.389262,-0.376705,-0.364148,-0.351591,-0.339035,-0.326478,-0.313921,-0.301364,-0.288807,-0.27625,-0.263694,-0.251137,-0.23858,-0.226023,-0.213466,-0.200909,-0.188353,-0.175796,-0.163239,-0.150682,-0.138125,-0.125568,-0.113012,-0.100455,-0.0878979,-0.075341,-0.0627842,-0.0502273,-0.0376705,-0.0251137,-0.0125568]






    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")









    



INFO: Saved animation to /Users/user/Programming/SIS/manual_burgers.gif






    Out[8]:



In [9]:

    
solution









    Out[9]:





200-element Array{Float64,1}:
  2.0      
  2.0      
  2.0      
  2.05479  
  2.07413  
  2.09334  
  2.11243  
  2.13138  
  2.15021  
  2.1689   
  2.18747  
  2.20592  
  2.22423  
  ⋮        
 -0.290221 
 -0.266061 
 -0.241895 
 -0.217723 
 -0.193546 
 -0.169364 
 -0.145177 
 -0.120986 
 -0.0967926
  0.0      
  0.0      
  0.0



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





DiffEqBase.ODESolution{Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,#rp},OrdinaryDiffEq.ImplicitEuler{0,true,Base.LinAlg.#lufact!},OrdinaryDiffEq.InterpolationData{#rp,Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.ImplicitEulerCache{Array{Float64,1},Array{Float64,1},Array{Float64,1},OrdinaryDiffEq.DiffCache{Float64,ForwardDiff.Dual{10,Float64}},Array{Float64,1},OrdinaryDiffEq.RHS_IE{#rp,Array{Float64,1},Float64,OrdinaryDiffEq.DiffCache{Float64,ForwardDiff.Dual{10,Float64}},Tuple{Int64},Base.OneTo{Int64}},NLsolve.DifferentiableMultivariateFunction,10}}}(Array{Float64,1}[[2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0  …  0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[2.0,1.99999,1.99999,1.99999,1.99998,1.99998,1.99998,1.99997,1.99997,1.99997  …  0.0118298,0.0104106,0.00906843,0.00779337,0.00657598,0.00540726,0.00427856,0.00318152,0.00210802,0.00105012],[1.99998,1.99997,1.99995,1.99993,1.99991,1.99989,1.99988,1.99986,1.99983,1.99981  …  0.0313001,0.0277692,0.0243667,0.0210773,0.017886,0.0147785,0.0117402,0.00875709,0.00581522,0.00290077],[1.99995,1.9999,1.99984,1.99979,1.99973,1.99968,1.99962,1.99956,1.9995,1.99943  …  0.0524691,0.0467794,0.0412299,0.0358049,0.0304886,0.0252654,0.0201196,0.0150356,0.009998,0.00499128],[1.99988,1.99977,1.99965,1.99953,1.9994,1.99928,1.99915,1.99901,1.99888,1.99873  …  0.0720364,0.0644206,0.0569346,0.0495643,0.0422951,0.035113,0.0280033,0.0209516,0.0139435,0.0069645],[1.99978,1.99956,1.99934,1.99912,1.99889,1.99866,1.99842,1.99818,1.99792,1.99767  …  0.0887823,0.079551,0.0704306,0.061409,0.0524743,0.0436141,0.0348163,0.0260684,0.0173581,0.00867278],[1.99965,1.99929,1.99893,1.99856,1.99819,1.99782,1.99744,1.99704,1.99664,1.99623  …  0.102548,0.0920042,0.0815507,0.0711782,0.0608768,0.0506367,0.0404477,0.0302998,0.0201827,0.0100862],[1.99947,1.99895,1.99841,1.99788,1.99733,1.99678,1.99622,1.99564,1.99506,1.99446  …  0.113611,0.102019,0.0904995,0.0790444,0.067646,0.0562965,0.0449879,0.0337122,0.022461,0.0112264],[1.99928,1.99855,1.99782,1.99708,1.99633,1.99557,1.9948,1.99402,1.99322,1.99241  …  0.12238,0.109961,0.0975981,0.0852864,0.0730191,0.0607903,0.0485935,0.0364225,0.0242709,0.0121322],[1.99906,1.99811,1.99715,1.99619,1.99522,1.99424,1.99324,1.99223,1.9912,1.99015  …  0.129261,0.116195,0.103173,0.0901895,0.0772407,0.0643216,0.0514273,0.038553,0.0256937,0.0128443]  …  [1.99246,1.98489,1.9773,1.96968,1.96204,1.95438,1.94669,1.93898,1.93124,1.92348  …  0.115819,0.104246,0.0926708,0.0810927,0.0695123,0.0579299,0.0463459,0.0347606,0.0231743,0.0115873],[1.99244,1.98486,1.97726,1.96963,1.96197,1.9543,1.9466,1.93887,1.93112,1.92335  …  0.115665,0.104108,0.0925478,0.080985,0.0694199,0.0578529,0.0462843,0.0347144,0.0231434,0.0115719],[1.99243,1.98484,1.97722,1.96958,1.96191,1.95422,1.94651,1.93877,1.93101,1.92322  …  0.115517,0.103974,0.0924286,0.0808806,0.0693304,0.0577783,0.0462246,0.0346695,0.0231136,0.0115569],[1.99242,1.98481,1.97718,1.96953,1.96185,1.95415,1.94642,1.93867,1.9309,1.9231  …  0.115372,0.103844,0.0923131,0.0807795,0.0692437,0.057706,0.0461667,0.0346261,0.0230846,0.0115425],[1.99241,1.98479,1.97715,1.96948,1.96179,1.95408,1.94634,1.93858,1.93079,1.92298  …  0.115233,0.103718,0.0922012,0.0806815,0.0691597,0.0576359,0.0461106,0.034584,0.0230566,0.0115284],[1.99239,1.98477,1.97711,1.96944,1.96173,1.95401,1.94626,1.93848,1.93069,1.92286  …  0.115097,0.103597,0.0920928,0.0805866,0.0690782,0.057568,0.0460563,0.0345433,0.0230294,0.0115148],[1.99238,1.98474,1.97708,1.96939,1.96168,1.95394,1.94618,1.9384,1.93059,1.92275  …  0.114966,0.103478,0.0919877,0.0804946,0.0689993,0.0575022,0.0460036,0.0345038,0.023003,0.0115017],[1.99237,1.98472,1.97705,1.96935,1.96162,1.95388,1.94611,1.93831,1.93049,1.92264  …  0.114839,0.103364,0.0918859,0.0804055,0.0689229,0.0574385,0.0459526,0.0344655,0.0229775,0.0114889],[1.99236,1.9847,1.97702,1.96931,1.96157,1.95381,1.94603,1.93823,1.93039,1.92254  …  0.114716,0.103253,0.0917873,0.0803191,0.0688488,0.0573768,0.0459032,0.0344284,0.0229528,0.0114765],[1.99235,1.98468,1.97699,1.96927,1.96152,1.95375,1.94596,1.93814,1.9303,1.92244  …  0.114597,0.103146,0.0916917,0.0802354,0.0687771,0.0573169,0.0458553,0.0343925,0.0229288,0.0114646]],[0.0,0.01,0.02,0.03,0.04,0.05,0.06,0.07,0.08,0.09  …  0.91,0.92,0.93,0.94,0.95,0.96,0.97,0.98,0.99,1.0],Array{Array{Float64,1},1}[Array{Float64,1}[[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0  …  0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0]],Array{Float64,1}[[0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0  …  0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],[-0.000301872,-0.000606932,-0.000917452,-0.00123575,-0.00156421,-0.0019053,-0.00226159,-0.00263577,-0.00303067,-0.00344928  …  1.18298,1.04106,0.906843,0.779337,0.657598,0.540726,0.427856,0.318152,0.210802,0.105012]],Array{Float64,1}[[-0.000301872,-0.000606932,-0.000917452,-0.00123575,-0.00156421,-0.0019053,-0.00226159,-0.00263577,-0.00303067,-0.00344928  …  1.18298,1.04106,0.906843,0.779337,0.657598,0.540726,0.427856,0.318152,0.210802,0.105012],[-0.00137847,-0.00276926,-0.00418046,-0.00562029,-0.00709715,-0.00861968,-0.0101968,-0.0118377,-0.013552,-0.0153496  …  1.94703,1.73586,1.52983,1.32839,1.13101,0.937121,0.746164,0.557558,0.37072,0.185065]],Array{Float64,1}[[-0.00137847,-0.00276926,-0.00418046,-0.00562029,-0.00709715,-0.00861968,-0.0101968,-0.0118377,-0.013552,-0.0153496  …  1.94703,1.73586,1.52983,1.32839,1.13101,0.937121,0.746164,0.557558,0.37072,0.185065],[-0.00348977,-0.00700618,-0.0105651,-0.0141825,-0.0178748,-0.0216585,-0.0255509,-0.0295694,-0.0337319,-0.0380571  …  2.1169,1.90101,1.68632,1.47277,1.26026,1.04869,0.837937,0.627851,0.418278,0.209052]],Array{Float64,1}[[-0.00348977,-0.00700618,-0.0105651,-0.0141825,-0.0178748,-0.0216585,-0.0255509,-0.0295694,-0.0337319,-0.0380571  …  2.1169,1.90101,1.68632,1.47277,1.26026,1.04869,0.837937,0.627851,0.418278,0.209052],[-0.0065179,-0.0130788,-0.0197053,-0.0264204,-0.0332474,-0.0402098,-0.0473317,-0.0546374,-0.0621519,-0.0699008  …  1.95673,1.76412,1.57047,1.37593,1.18065,0.984757,0.78837,0.591601,0.394553,0.197322]],Array{Float64,1}[[-0.0065179,-0.0130788,-0.0197053,-0.0264204,-0.0332474,-0.0402098,-0.0473317,-0.0546374,-0.0621519,-0.0699008  …  1.95673,1.76412,1.57047,1.37593,1.18065,0.984757,0.78837,0.591601,0.394553,0.197322],[-0.0100771,-0.0202122,-0.0304321,-0.0407637,-0.0512341,-0.0618708,-0.0727015,-0.0837543,-0.0950578,-0.106641  …  1.6746,1.51304,1.34959,1.18448,1.01791,0.850116,0.681304,0.511682,0.341455,0.170828]],Array{Float64,1}[[-0.0100771,-0.0202122,-0.0304321,-0.0407637,-0.0512341,-0.0618708,-0.0727015,-0.0837543,-0.0950578,-0.106641  …  1.6746,1.51304,1.34959,1.18448,1.01791,0.850116,0.681304,0.511682,0.341455,0.170828],[-0.013713,-0.0274962,-0.0413767,-0.0553819,-0.0695396,-0.0838774,-0.0984234,-0.113206,-0.128253,-0.143593  …  1.3766,1.24532,1.11201,0.976918,0.840258,0.702258,0.563142,0.423136,0.282461,0.141342]],Array{Float64,1}[[-0.013713,-0.0274962,-0.0413767,-0.0553819,-0.0695396,-0.0838774,-0.0984234,-0.113206,-0.128253,-0.143593  …  1.3766,1.24532,1.11201,0.976918,0.840258,0.702258,0.563142,0.423136,0.282461,0.141342],[-0.0170495,-0.0341773,-0.0514085,-0.0687683,-0.0862818,-0.103974,-0.121871,-0.139999,-0.158381,-0.177045  …  1.10631,1.00151,0.894885,0.786618,0.676916,0.565982,0.45402,0.341236,0.227834,0.11402]],Array{Float64,1}[[-0.0170495,-0.0341773,-0.0514085,-0.0687683,-0.0862818,-0.103974,-0.121871,-0.139999,-0.158381,-0.177045  …  1.10631,1.00151,0.894885,0.786618,0.676916,0.565982,0.45402,0.341236,0.227834,0.11402],[-0.0198466,-0.0397764,-0.0598102,-0.0799693,-0.100275,-0.120748,-0.141409,-0.16228,-0.183381,-0.204734  …  0.876823,0.794123,0.709859,0.624199,0.537315,0.449378,0.36056,0.271038,0.180987,0.0905818]],Array{Float64,1}[[-0.0198466,-0.0397764,-0.0598102,-0.0799693,-0.100275,-0.120748,-0.141409,-0.16228,-0.183381,-0.204734  …  0.876823,0.794123,0.709859,0.624199,0.537315,0.449378,0.36056,0.271038,0.180987,0.0905818],[-0.0219971,-0.0440793,-0.0662628,-0.0885637,-0.110998,-0.133583,-0.156333,-0.179265,-0.202395,-0.225737  …  0.688168,0.623458,0.557461,0.490315,0.422158,0.353132,0.28338,0.213046,0.142274,0.0712102]]  …  Array{Float64,1}[[-0.0013798,-0.00276357,-0.00415098,-0.00554167,-0.00693528,-0.00833147,-0.00972986,-0.0111301,-0.0125319,-0.0139347  …  -0.0163516,-0.0147292,-0.0131028,-0.0114728,-0.00983967,-0.00820384,-0.00656576,-0.0049259,-0.00328468,-0.00164256],[-0.00133694,-0.00267774,-0.00402206,-0.00536955,-0.00671988,-0.0080727,-0.00942767,-0.0107844,-0.0121426,-0.0135019  …  -0.015844,-0.0142719,-0.012696,-0.0111166,-0.00953418,-0.00794914,-0.00636192,-0.00477297,-0.0031827,-0.00159157]],Array{Float64,1}[[-0.00133694,-0.00267774,-0.00402206,-0.00536955,-0.00671988,-0.0080727,-0.00942767,-0.0107844,-0.0121426,-0.0135019  …  -0.015844,-0.0142719,-0.012696,-0.0111166,-0.00953418,-0.00794914,-0.00636192,-0.00477297,-0.0031827,-0.00159157],[-0.00129543,-0.00259459,-0.00389716,-0.00520281,-0.00651121,-0.00782202,-0.00913491,-0.0104495,-0.0117656,-0.0130827  …  -0.015352,-0.0138288,-0.0123018,-0.0107714,-0.00923813,-0.00770231,-0.00616438,-0.00462476,-0.00308387,-0.00154215]],Array{Float64,1}[[-0.00129543,-0.00259459,-0.00389716,-0.00520281,-0.00651121,-0.00782202,-0.00913491,-0.0104495,-0.0117656,-0.0130827  …  -0.015352,-0.0138288,-0.0123018,-0.0107714,-0.00923813,-0.00770231,-0.00616438,-0.00462476,-0.00308387,-0.00154215],[-0.00125521,-0.00251403,-0.00377617,-0.00504128,-0.00630906,-0.00757918,-0.0088513,-0.0101251,-0.0114003,-0.0126765  …  -0.0148752,-0.0133993,-0.0119197,-0.0104369,-0.00895122,-0.0074631,-0.00597293,-0.00448113,-0.0029881,-0.00149425]],Array{Float64,1}[[-0.00125521,-0.00251403,-0.00377617,-0.00504128,-0.00630906,-0.00757918,-0.0088513,-0.0101251,-0.0114003,-0.0126765  …  -0.0148752,-0.0133993,-0.0119197,-0.0104369,-0.00895122,-0.0074631,-0.00597293,-0.00448113,-0.0029881,-0.00149425],[-0.00121624,-0.00243599,-0.00365895,-0.00488479,-0.00611321,-0.0073439,-0.00857654,-0.00981082,-0.0110464,-0.012283  …  -0.0144132,-0.0129831,-0.0115495,-0.0101127,-0.00867319,-0.00723128,-0.0057874,-0.00434194,-0.00289529,-0.00144784]],Array{Float64,1}[[-0.00121624,-0.00243599,-0.00365895,-0.00488479,-0.00611321,-0.0073439,-0.00857654,-0.00981082,-0.0110464,-0.012283  …  -0.0144132,-0.0129831,-0.0115495,-0.0101127,-0.00867319,-0.00723128,-0.0057874,-0.00434194,-0.00289529,-0.00144784],[-0.00117849,-0.00236039,-0.00354538,-0.00473318,-0.00592348,-0.00711597,-0.00831035,-0.00950631,-0.0107035,-0.0119017  …  -0.0139654,-0.0125798,-0.0111907,-0.00979857,-0.00840376,-0.00700664,-0.00560762,-0.00420706,-0.00280534,-0.00140286]],Array{Float64,1}[[-0.00117849,-0.00236039,-0.00354538,-0.00473318,-0.00592348,-0.00711597,-0.00831035,-0.00950631,-0.0107035,-0.0119017  …  -0.0139654,-0.0125798,-0.0111907,-0.00979857,-0.00840376,-0.00700664,-0.00560762,-0.00420706,-0.00280534,-0.00140286],[-0.00114192,-0.00228714,-0.00343536,-0.00458629,-0.00573965,-0.00689513,-0.00805245,-0.00921129,-0.0103714,-0.0115324  …  -0.0135315,-0.0121889,-0.010843,-0.00949414,-0.00814266,-0.00678896,-0.0054334,-0.00407635,-0.00271819,-0.00135928]],Array{Float64,1}[[-0.00114192,-0.00228714,-0.00343536,-0.00458629,-0.00573965,-0.00689513,-0.00805245,-0.00921129,-0.0103714,-0.0115324  …  -0.0135315,-0.0121889,-0.010843,-0.00949414,-0.00814266,-0.00678896,-0.0054334,-0.00407635,-0.00271819,-0.00135928],[-0.00110649,-0.00221616,-0.00332876,-0.00444398,-0.00556154,-0.00668117,-0.00780258,-0.00892546,-0.0100495,-0.0111745  …  -0.0131111,-0.0118102,-0.0105061,-0.00919915,-0.00788966,-0.00657802,-0.00526457,-0.00394969,-0.00263373,-0.00131704]],Array{Float64,1}[[-0.00110649,-0.00221616,-0.00332876,-0.00444398,-0.00556154,-0.00668117,-0.00780258,-0.00892546,-0.0100495,-0.0111745  …  -0.0131111,-0.0118102,-0.0105061,-0.00919915,-0.00788966,-0.00657802,-0.00526457,-0.00394969,-0.00263373,-0.00131704],[-0.00107215,-0.0021474,-0.00322547,-0.00430609,-0.00538898,-0.00647387,-0.00756048,-0.00864853,-0.00973774,-0.0108278  …  -0.0127037,-0.0114432,-0.0101797,-0.00891329,-0.0076445,-0.00637361,-0.00510098,-0.00382696,-0.00255189,-0.00127612]],Array{Float64,1}[[-0.00107215,-0.0021474,-0.00322547,-0.00430609,-0.00538898,-0.00647387,-0.00756048,-0.00864853,-0.00973774,-0.0108278  …  -0.0127037,-0.0114432,-0.0101797,-0.00891329,-0.0076445,-0.00637361,-0.00510098,-0.00382696,-0.00255189,-0.00127612],[-0.00103889,-0.00208078,-0.0031254,-0.0041725,-0.00522179,-0.00627302,-0.00732592,-0.00838021,-0.00943562,-0.0104919  …  -0.0123089,-0.0110876,-0.00986331,-0.00863631,-0.00740694,-0.00617555,-0.00494246,-0.00370803,-0.00247258,-0.00123646]],Array{Float64,1}[[-0.00103889,-0.00208078,-0.0031254,-0.0041725,-0.00522179,-0.00627302,-0.00732592,-0.00838021,-0.00943562,-0.0104919  …  -0.0123089,-0.0110876,-0.00986331,-0.00863631,-0.00740694,-0.00617555,-0.00494246,-0.00370803,-0.00247258,-0.00123646],[-0.00100666,-0.00201623,-0.00302845,-0.00404305,-0.0050598,-0.00607842,-0.00709865,-0.00812023,-0.0091429,-0.0101664  …  -0.0119264,-0.010743,-0.00955678,-0.00836791,-0.00717675,-0.00598362,-0.00478886,-0.00359279,-0.00239574,-0.00119803]]],DiffEqBase.ODEProblem{Array{Float64,1},Float64,true,#rp}(rp,[2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0,2.0  …  0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0,0.0],(0.0,1.0)),OrdinaryDiffEq.ImplicitEuler{0,true,Base.LinAlg.#lufact!}(lufact!),OrdinaryDiffEq.InterpolationData,true,0)



In [11]:

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









    



(Base.OneTo(101),)






    Out[11]:



In [12]:

    
sol[end,:]









    Out[12]:





200-element Array{Float64,1}:
 1.99235  
 1.98468  
 1.97699  
 1.96927  
 1.96152  
 1.95375  
 1.94596  
 1.93814  
 1.9303   
 1.92244  
 1.91455  
 1.90663  
 1.89869  
 ⋮        
 0.13749  
 0.126045 
 0.114597 
 0.103146 
 0.0916917
 0.0802354
 0.0687771
 0.0573169
 0.0458553
 0.0343925
 0.0229288
 0.0114646



In [ ]: