## Solving *steady* state series equations
\begin{equation}
\Psi _{\eta \eta \tau }^{(n)}=\frac{1}{R_s}\Psi_{\eta \eta \eta \eta
}^{(n)}+ \sum\limits_{i = 0}^n{(2i+1)\left[ \Psi ^{(i)} \Psi _{\eta \eta
\eta }^{(n-i)} -\Psi_{\eta}^{(n-i)} \Psi _{\eta \eta}^{(i)}\right] } \\
\quad \Psi^{(n)}(\eta=\mp 1) = 0, \qquad \Psi^{(n)}_\eta(\eta=\mp 1) = \frac{(-1)^{n+1}}{(2n+1)!}, \qquad
-1 \le \eta \le1
\end{equation}

The steady state equation above will be solved using the two point boundary value solver of Julia and shooting methods.


```julia
using DifferentialEquations
using Plots
using Sundials
using BenchmarkTools
using Calculus
using DelimitedFiles
using Interpolations
using LaTeXStrings
using FFTW
```

## Solving the leading order equation $\Psi(\eta)^{(0)}$ steady state solution

### Low Reynolds number $R<17$
For small values of the Reynolds number ($R<17)$, the first order approximation $\Psi^{(0)}$ can be calculated fairly quickly using standard boundary value solvers (memory allocation is around 15GB at R=50 using MIRK4 two point boundary value solver). With a good initial guess value and correct values for ```abstol``` and ```reltol``` the Shooting method also works very well with all solvers e.g. Vern7()

### Intermediate Reynolds numbers $17<R<65$
Despite the equations being unstable to computation for $R>65$, they are still solvable for $17<R<65$ and the only reliable way to solve for $\Psi^{(0)}$ is to use the ```TwoPointBVProblem``` or Shooting with GeneralMirk4/Rosenbrock23. The latter solver requirement is due to the equations becoming stiff.


```julia
sol_psi = Array{ODESolution,1}(undef, 20) #Declaring an array for every \psi(\eta)^{(n)} solution
series_term = 0
Reynolds_number = 1.2345

p=(Reynolds_number, series_term)
const etaspan=(-1.0,1.0)
ψ0= [0.0, -1.0, 3, -3] #This contains psi^{(0)}, psi^{(0)}_eta, psi^{(0)}_eta_eta, psi^{(0)}_eta_eta_eta


"""This is the differential equation for the ``\\psi(\\eta)^{(0)}`` term in the series.
``\\Psi_{\\eta \\eta \\eta \\eta
}^{(n)} = R_s( \\Psi_{\\eta}^{(0)} \\Psi _{\\eta \\eta}^{(0)} - \\Psi ^{(0)} \\Psi _{\\eta \\eta
\\eta }^{(0)})``
"""
function steady_diffeq_psi0!(dψ, ψ, p, η)
  dψ[1] = ψ[2]
  dψ[2] = ψ[3]
  dψ[3] = ψ[4]
  dψ[4] = (ψ[2]*ψ[3]-ψ[4]*ψ[1])*p[1]
end

"""This is the differential equation for the ``\\psi(\\eta)^{(n)}`` term in the series.
``\\Psi_{\\eta \\eta \\eta \\eta
}^{(n)} = R_s\\sum\\limits_{i = 0}^n{(2i+1)( \\Psi_{\\eta}^{(n-i)} \\Psi _{\\eta \\eta}^{(i)} - \\Psi ^{(i)} \\Psi _{\\eta \\eta
\\eta }^{(n-i)})}``
"""
function steady_diffeq_series!(dψ, ψ, p, η)
    solut_psi0 = p[3]::ODESolution #This avoids the use of Global sol_psi0 and should speed up the code, this variable type checked
    solut_psi  = p[4]::Vector{ODESolution}  #This avoids the use of Global sol_psi and should speed up the code,  this variable type checked
    sum_of_nonlinears = 0.0
    n = round(Int, p[2]) # This rounding off prevents errors are thrown when p vector is changed by Julia from Int to Floats
    for i= 1:(n-1)
        sum_of_nonlinears += (2*i+1)*(solut_psi[n-i](η)[2]*solut_psi[i](η)[3]
            - solut_psi[i](η)[1]*solut_psi[n-i](η)[4])
    end
    dψ[1] = ψ[2]::Float64
    dψ[2] = ψ[3]::Float64
    dψ[3] = ψ[4]::Float64
    dψ[4] = (ψ[2]*solut_psi0(η)[3]-solut_psi0(η)[1]*ψ[4])*p[1] +
           (2*p[2]+1)*(solut_psi0(η)[2]*ψ[3]-ψ[1]*solut_psi0(η)[4])*p[1] +
            sum_of_nonlinears*p[1]
end

"""
This is the Jacobian function for the steady ``\\psi(\\eta)^{(0)}`` differential equation
"""
function jacobian_steadypsi0(J,ψ,p,η)

    J[1,1] = J[1,3] = J[1,4] =  0.0;
    J[1,2] = 1.0
  
    J[2,1] = J[2,2] = J[2,4] =  0.0;
    J[2,3] = 1.0  
    
    J[3,1] = J[3,2] = J[3,3] =  0.0;  
    J[3,4] = 1.0  
    
    J[3,1] = -ψ[4]*p[1]
    J[3,2] =  ψ[3]*p[1]
    J[3,3] =  ψ[2]*p[1]
    J[3,4] = -ψ[1]*p[1]    
  nothing
end


"""
This is the boundary condition for the nth term  ``\\psi(\\eta)^{(n)}``. The p-vector supplies the correct value 
of n in the 2nd position i.e. p[2] = n
"""
function bc!(residual, ψ, p, η)  # psi[1] is the beginning of the etaspan, and psi[end] is the ending
    residual[1] = ψ[1][1]         # The psi[1] (i.e,. psi^{0}) solution at the beginning of the time span should be 0
    residual[2] = ψ[1][2] - (-1)^(p[2]+1)/factorial(2*p[2]+1)       #First derivative of  should be -1 at first time step
    residual[3] = ψ[end][1]       # the solution at the end of the time span should be 0
    residual[4] = ψ[end][2] - (-1)^(p[2]+1)/factorial(2*p[2]+1)      #First derivative should be -1 at end time step
end

function psi_smallr(i, η ,r) 
    return((-1.0)^(i+1)*η*(η^2-1)/(2.0*factorial(2*i+1)) + r*(-1)^i*(2.0)^(2*i-3)*η*(η^2-1)^2.0*(2.0+η^2)/(35.0*factorial(2*i+1))+r^2.0*η*(η^2-1)^2*((-591 -2294*η^2 + 161*η^4 + 1428*η^6)*(-1.0)^i + 
                        3.0*(423 + 166*η^2 - 553*η^4 - 84*η^6)*(-1)^i*(3)^(2*i+1))/(62092800*factorial(2*i+1)))
end

function psi_smallr_d1(i, η,  r ) 
    return((-1.0)^(i+1)*(3.0*η^2-1)/(2.0*factorial(2*i+1)) + r*(-1)^i*(2.0)^(2*i-3)*(7*η^6-9*η^2+2)/(35.0*factorial(2*i+1)))
end

function psi_smallr_d2(i, η, r) 
    return((-1.0)^(i+1)*(6*η)/(2.0*factorial(2*i+1)) + r*(-1)^i*(2.0)^(2*i-3)*(42*η^5-18*η)/(35.0*factorial(2*i+1))+
    r^2.0*(3*(-197 - 1112*η^2 + 6930*η^4 - 2772*η^6 - 8085*η^8 + 5236*η^10)*(-1.0)^i +
            (1269 - 6120*η^2 - 6930*η^4 + 24948*η^6 - 10395*η^8 - 2772*η^10)*(-1)^i*(3)^(2*i+1)) /(62092800*factorial(2*i+1)))
end

function psi_smallr_d3(i, η,  r) 
    return((-1.0)^(i+1)*(6)/(2.0*factorial(2*i+1)) + r*(-1)^i*(2.0)^(2*i-3)*(210*η^4-18)/(35.0*factorial(2*i+1)))
end

function psi_bigr(i, η ,r) 
    return( (-1)^i *sin(pi*η)/(pi*factorial(2*i+1)))
end

function psi_bigr_d1(i, η ,r) 
    return( (-1)^i *cos(pi*η)/factorial(2*i+1))
end

function psi_bigr_d3(i, η ,r) 
    return( -(-1)^i*pi^2*cos(pi*η)/factorial(2*i+1))
end

psi0_d2_stored = [2.91, 2.57,  2.18, 1.82,  (1.5-0.3),       1.02,  0.847,   0.71, 0.599, 0.511, 0.442,  0.38,  0.33]
psi0_d3_stored = [-2.32,0.22,  2.98, 5.2,  (7.02+5.0), (9.29+2.7), 9.9598,  (10.42+1), 10.72, 10.93, 11.05, 11.14, 11.19]
Reynolds_values_stored = [1.0,5,10,  15,           20,         30,     35,     40,    45,    50,    55,    60,    65]
itp_d2 = Interpolations.interpolate((Reynolds_values_stored,), psi0_d2_stored , Gridded(Linear()))
itp_d3 = Interpolations.interpolate((Reynolds_values_stored,), psi0_d3_stored, Gridded(Linear()))

function psi0_bc_initialguess(Reynolds_number)
    if (Reynolds_number <=65)
        return[0.0, -1.0, itp_d2(Reynolds_number), itp_d3(Reynolds_number)] #This contains psi, psi_eta, psi_eta_eta, psi_eta_eta_eta
    else
        #return([0.0, -1.0, 0.33, 12])
         return[0.0, psi_bigr_d1(0, 0.0 ,Reynolds_number), 0.0, psi_bigr_d3(0, 0.0 ,Reynolds_number)] #This contains ψ, ψ_η, ψ_ηη, ψ_ηηη
    end
end

#The code below is for testing purposes only, it treats the psi problem as an IVP starting at eta=-1
#p=(Reynolds_number, 0)
#ψ0= [0.0, -1.0, 0.02, 2.9836060306009644] 
#prob = ODEProblem(steady_diffeq_psi0!,ψ0,etaspan,p , jac = jacobian_steadypsi0)
#sol =  solve(prob,reltol=1e-6)
#@benchmark solve(prob,reltol=1e-6,save_everystep=true)
```




    psi0_bc_initialguess (generic function with 1 method)




```julia
Reynolds_number = 10000.0
series_term = 0
p=(Reynolds_number, series_term)
ψ0 = psi0_bc_initialguess(Reynolds_number)
bvp_psi0_2point = TwoPointBVProblem(steady_diffeq_psi0!, bc!, ψ0, etaspan, p , jac = jacobian_steadypsi0)
sol_psi0 =  @time solve(bvp_psi0_2point, alg_hints = [:stiff], MIRK4(),dt=0.01) #Very accurate solver

```

     15.856637 seconds (220.40 M allocations: 20.656 GiB, 16.86% gc time, 0.20% compilation time)





    retcode: Success
    Interpolation: 1st order linear
    t: 201-element Vector{Float64}:
     -1.0
     -0.99
     -0.98
     -0.97
     -0.96
     -0.95
     -0.94
     -0.93
     -0.92
     -0.91
     -0.9
     -0.89
     -0.88
      ⋮
      0.89
      0.9
      0.91
      0.92
      0.93
      0.94
      0.95
      0.96
      0.97
      0.98
      0.99
      1.0
    u: 201-element Vector{Vector{Float64}}:
     [-5.438454421878129e-37, -1.0, 0.00012300546589258518, 9.906643500646569]
     [-0.009998343261789869, -0.9995036509608582, 0.09912384398995155, 9.892868477000627]
     [-0.019986775490410332, -0.998018068540089, 0.19795957798109, 9.873031167561809]
     [-0.02995541373759932, -0.9955452515082693, 0.29655721166609045, 9.845061137444274]
     [-0.039894398957369744, -0.992088002421085, 0.39483122651067376, 9.808233581305148]
     [-0.049793904502366176, -0.9876500063959832, 0.49269133580438895, 9.762247958064613]
     [-0.059644145131804965, -0.9822358630244032, 0.5900452162596184, 9.706974606811162]
     [-0.06943538623439384, -0.9758511000072576, 0.6867997078941593, 9.642367041490468]
     [-0.07915795313930095, -0.9685021780843683, 0.7828614416842025, 9.568426041749166]
     [-0.08880224044997845, -0.9601964910890403, 0.8781372189316826, 9.485182971882717]
     [-0.09835872136191663, -0.9509423628802954, 0.9725342699164605, 9.392691231300525]
     [-0.1078179569380116, -0.9407490420413301, 1.0659604483402518, 9.291021509761542]
     [-0.11717060532189127, -0.9296266948334767, 1.1583243889991424, 9.180258971108556]
     ⋮
     [0.10781795693803839, -0.9407490420420476, -1.065960448327419, 9.291021509653053]
     [0.09835872136193688, -0.9509423628808888, -0.9725342699047156, 9.392691231189698]
     [0.08880224044999327, -0.9601964910895239, -0.8781372189210467, 9.485182971771357]
     [0.0791579531393115, -0.9685021780847505, -0.7828614416746874, 9.568426041636064]
     [0.06943538623440214, -0.9758511000075507, -0.6867997078857929, 9.642367041375119]
     [0.05964414513181105, -0.9822358630246191, -0.5900452162524249, 9.706974606693553]
     [0.04979390450237011, -0.9876500063961322, -0.4926913357983836, 9.762247957946519]
     [0.0398943989573728, -0.9920880024211807, -0.3948312265058404, 9.80823358118756]
     [0.029955413737601678, -0.995545251508323, -0.29655721166244514, 9.8450611373168]
     [0.019986775490412237, -0.9980180685401135, -0.197959577978685, 9.873031167453409]
     [0.009998343261791635, -0.9995036509608635, -0.09912384398864327, 9.892868476887353]
     [-7.239780535325076e-36, -1.0, -0.0001230054660869085, 9.906643500568554]




```julia
if(p[1] < 52.5)
    run(`cc -ansi steady_d4_julia.c -o steadyd4julia -lm `)
    io = open("Reynolds_number.txt", "w");
    write(io, string(p[1]));
    close(io);
    sleep(1)
    run(`./steadyd4julia $p\[1\]`)
    MatrixPsis=readdlm("terms_c.data",Float64,);
end
```


```julia
print(varinfo(r"sol_psi0"))
pyplot()
plot(sol_psi0,vars=(1), linecolor=:red, label = "Julia solution",xlabel=L"\eta", ylabel=L"\Psi^{(0)}")
if(p[1] < 52.5)
    plot!(MatrixPsis[:,1], MatrixPsis[:, 2], linestyle=:dot, linecolor=:green, label = "C code solution")
    plot!( x-> psi_smallr(0, x, p[1]), -1,1,linestyle=:dot, linecolor=:orange, label = "Small R analytical",title="R = "*string(p[1])*L",   \Psi^{(0)}", legend=:bottomright)
end

plot!( x-> psi_bigr(0, x, p[1]), -1,1,linestyle=:dash, linecolor=:blue, label = "Large R analytical",title="R = "*string(p[1])*L",   \Psi^{(0)}", legend=:bottomright, yguidefontrotation=-90)

```


```julia
psi_plot = Array{Plots.Plot{Plots.PyPlotBackend},1}()
u0 =[]
max_terms = 5
for i = 1:max_terms
    series_term = i
    p=(Reynolds_number, series_term, sol_psi0, sol_psi)
    if (Reynolds_number > 20 )
        push!(u0, [0.0, psi_smallr_d1(series_term, -1.0, p[1]), 0,  pi^2*psi_smallr_d1(series_term, -1.0, p[1])])  
    else
        push!(u0, [0.0, psi_smallr_d1(series_term, -1.0, p[1]), psi_smallr_d2(series_term, -1.0, p[1]), psi_smallr_d3(series_term, -1.0, p[1])])  
    end
    
    bvp_sol= TwoPointBVProblem(steady_diffeq_series!, bc!, u0[i], etaspan, p)
    sol_psi[ i ] = @time solve(bvp_sol, alg_hints=[:stiff], MIRK4(), dt=0.04)
    plot(sol_psi[i],vars=(1), linecolor=:red, label = "Julia solution",xlabel=L"\eta", ylabel=L"\Psi^{(%$i)}", size=(800,400))
    if(p[1] < 35)
        plot!(MatrixPsis[:,1],MatrixPsis[:,i+2], linecolor=:green, label = "C code solution")
    end
    if(p[1] < 35  )
            plot!(x-> psi_smallr(i, x, p[1]),-1,1,linestyle=:dot, linecolor=:orange, label = "Small R analytical",title="R = "*string(p[1])*L",   \Psi^{(%$i)}", legend=:outertopright)
    end
    push!(psi_plot,plot!(x-> psi_bigr(i, x, p[1]),-1,1,linestyle=:dash, linecolor=:blue, label = "Large R analytical",title="R = "*string(p[1])*L",   \Psi^{(%$i)}", legend=:outertopright, yguidefontrotation=-90))
end
```

      0.489074 seconds (7.21 M allocations: 331.166 MiB, 10.30% gc time, 14.26% compilation time)
      0.673325 seconds (12.80 M allocations: 601.316 MiB, 13.20% gc time)
      0.974780 seconds (19.91 M allocations: 942.041 MiB, 11.86% gc time)
      1.280404 seconds (27.03 M allocations: 1.253 GiB, 12.02% gc time)
      4.370773 seconds (93.95 M allocations: 4.363 GiB, 12.27% gc time)



```julia
for i in 1:max_terms
          display(psi_plot[i])
          sleep(1)
 end

```

## Exploring with Fourier Transforms


```julia
alist= [sol_psi0.u[i][1] for i =1:length(sol_psi0.u)]
x = range(-1,stop=1,length=length(sol_psi[5].u))
clist= [sol_psi0.u[i][3] for i =1:length(sol_psi0.u)]


clist[25]/alist[25]
```




    -11.30626039872988




```julia
N=length(sol_psi[5].u)
Fy = fft(alist)[1:N÷2]
ak =  2/N * real.(Fy)
bk = -2/N * imag.(Fy)  # fft sign convention
ak[1] = ak[1]/2;
```


```julia
bk*10^5
```




    25-element Vector{Float64}:
     -0.0
      5.588135150784316
     -1.4677891682918485
     -1.7602036153999576
     -0.8297559692080744
     -0.016353052152772343
      0.25189645884487194
      0.21082570799079534
      0.11893810212564905
      0.06150421055572241
      0.03369089940397449
      0.02000129490045747
      0.012565369887931823
      0.008159797258145592
      0.0053920674968345435
      0.0035852750426532366
      0.002376731460287785
      0.0015572299338542958
      0.0009991303702340953
      0.0006209769328401041
      0.0003687704513288578
      0.00020545369856001777
      0.00010470660852719787
      4.7114774678100995e-5
      1.7690957441206614e-5




```julia

```




    

    



## Solving the for $\Psi^{0}$ using shooting method


```julia
Reynolds_number =65.0
series_term = 0
ψ0= psi0_bc_initialguess(Reynolds_number)
p=(Reynolds_number, series_term)
bvp_shooting = BVProblem(steady_diffeq_psi0!, bc!, ψ0, etaspan,p , jac=jacobian_steadypsi0)
sol_shooting = @time solve(bvp_shooting, alg_hints=[:stiff], reltol=1e-5, abstol=1e-5 , GeneralMIRK4(),dt=0.01) #Excellent solver
#sol_shooting = @time solve(bvp_shooting, alg_hints=[:stiff], reltol=1e-8, abstol=1e-6 , Shooting(CVODE_BDF())) #Ok but need to have tolerances adjusted some errors on boundary
#sol_shooting = @time solve(bvp_shooting, alg_hints=[:stiff], reltol=1e-8, abstol=1e-6 , Shooting(Rodas4P())) #Ok but need to have tolerances adjusted some errors on boundary
#sol_shooting = @time solve(bvp_shooting, alg_hints=[:stiff], reltol=1e-9, abstol=1e-9 , Shooting(RadauIIA5())) #Good, FAST, STABLE but some errors on boundary
#sol_shooting = @time solve(bvp_shooting, alg_hints=[:stiff], reltol=1e-8, abstol=1e-6 , Shooting(TRBDF2()))
#sol_shooting = @time solve(bvp_shooting, alg_hints=[:stiff], reltol=1e-6, abstol=1e-6 , Shooting(Rosenbrock23()))


```

      5.704660 seconds (78.20 M allocations: 7.331 GiB, 18.45% gc time)





    retcode: Success
    Interpolation: 1st order linear
    t: 201-element Vector{Float64}:
     -1.0
     -0.99
     -0.98
     -0.97
     -0.96
     -0.95
     -0.94
     -0.93
     -0.92
     -0.91
     -0.9
     -0.89
     -0.88
      ⋮
      0.89
      0.9
      0.91
      0.92
      0.93
      0.94
      0.95
      0.96
      0.97
      0.98
      0.99
      1.0
    u: 201-element Vector{Vector{Float64}}:
     [0.0, -1.0, 0.3381169611051619, 11.193291745248118]
     [-0.009981237755586674, -0.9960628262890472, 0.4489518580668842, 10.973737934551135]
     [-0.0199175986148863, -0.991028279732876, 0.5575913604840704, 10.754082323461496]
     [-0.02979821866776514, -0.9849183318338913, 0.6640307658781293, 10.53361026174277]
     [-0.039612454067938174, -0.9777550347096055, 0.7682588601435602, 10.311730960405413]
     [-0.049349882144788734, -0.9695605803143844, 0.8702590678114968, 10.087960306035388]
     [-0.05900030305268773, -0.9603573489715499, 0.9700104441267903, 9.861906327821275]
     [-0.06855374185831681, -0.9501679486729377, 1.067488533169761, 9.633256872063251]
     [-0.07800045098008439, -0.9390152463814293, 1.1626661122064241, 9.401769117992206]
     [-0.08733091290498911, -0.9269223923874316, 1.2555138391191516, 9.167260632923332]
     [-0.0965358431181567, -0.9139128386157033, 1.3460008170328064, 8.929601717338844]
     [-0.10560619318841406, -0.900010351650049, 1.434095087958115, 8.688708833356118]
     [-0.11453315396047871, -0.8852390211341102, 1.519764065416849, 8.44453894522796]
     ⋮
     [0.1056061931885953, -0.9000103516534113, -1.4340950879367813, 8.68870883360754]
     [0.09653584311830656, -0.913912838618844, -1.3460008170090008, 8.929601717570556]
     [0.08733091290511165, -0.9269223923903188, -1.2555138390930447, 9.167260633133235]
     [0.0780004509801867, -0.939015246384057, -1.162666112178576, 9.401769118182791]
     [0.06855374185839476, -0.9501679486752653, -1.0674885331402812, 9.633256872233433]
     [0.059000303052743154, -0.9603573489735784, -0.9700104440955563, 9.861906327970678]
     [0.04934988214482538, -0.9695605803161094, -0.870259067778888, 10.087960306165272]
     [0.039612454067963106, -0.9777550347109999, -0.7682588601097937, 10.311730960512747]
     [0.0297982186677809, -0.9849183318349432, -0.6640307658433727, 10.533610261829677]
     [0.01991759861489311, -0.9910282797335839, -0.5575913604487551, 10.754082323524418]
     [0.009981237755591175, -0.996062826289411, -0.448951858031075, 10.973737934592188]
     [1.1692924532960299e-33, -1.0, -0.33811696106907635, 11.193291745265805]




```julia
plot(sol_shooting, vars=(1))
```


```julia

```

046759a30361724648710f3e9a4a971cb97fba7d
ipynb
Jupyter Notebook
Steadystreaming_series_solution_BVP.ipynb
gsagoo/SolvingDifferentialEquations
63c2b231b45ad64206a11824198b3d5fe20c062b
Unlicense
Steadystreaming_series_solution_BVP.ipynb
gsagoo/SolvingDifferentialEquations
63c2b231b45ad64206a11824198b3d5fe20c062b
Unlicense
Steadystreaming_series_solution_BVP.ipynb
gsagoo/SolvingDifferentialEquations
63c2b231b45ad64206a11824198b3d5fe20c062b
Unlicense
Qwen/Qwen-72B
1. YES
2. YES