Giving a try, and they are tons of mistake because I’m not 100% sure about the eq and bc,
julia> using ModelingToolkit,MethodOfLines,DomainSets
julia> @parameters k,ρ
2-element Vector{Num}:
k
ρ
julia> @variables V(..)
1-element Vector{Symbolics.CallWithMetadata{SymbolicUtils.FnType{Tuple, Real}, Base.ImmutableDict{DataType, Any}}}:
V⋆
julia> Dρ = Differential(ρ)
Differential(ρ)
julia> D2ρ = Dρ^2
Differential(ρ) ∘ Differential(ρ)
julia> Dk = Differential(k)
Differential(k)
julia> kmin = -4e5
-400000.0
julia> kmax= 3e5
300000.0
julia> rhomin = 0
0
julia> rhomax = 2e4
20000.0
julia> T= 1
1
julia> N = 2
2
julia> eq = [k*Dk(V(k,ρ)) ~ T*k^3*(k^2/(k^2+Dρ(V(k,ρ)) + 2*ρ*D2ρ(V(k,ρ))) + (N-1)*k^2/(k^2 + Dρ(V(k,ρ))))]
1-element Vector{Equation}:
k*Differential(k)(V(k, ρ)) ~ (k^3)*((k^2) / (Differential(ρ)(V(k, ρ)) + k^2 + 2Differential(ρ)(Differential(ρ)(V(k, ρ)))*ρ) + (k^2) / (Differential(ρ)(V(k, ρ)) + k^2))
julia> bcs = [V(kmin,ρ)~λ*ρ, V(k,0)~0,Dρ(V(k,0)) ~ λ ]
3-element Vector{Equation}:
V(-400000.0, ρ) ~ 2ρ
V(k, 0) ~ 0
Differential(ρ)(V(k, 0)) ~ 2
julia> @named pdesys = PDESystem(eq,bcs,domains,[k,ρ],[V(k,ρ)])
PDESystem
Equations: Equation[k*Differential(k)(V(k, ρ)) ~ (k^3)*((k^2) / (Differential(ρ)(V(k, ρ)) + k^2 + 2Differential(ρ)(Differential(ρ)(V(k, ρ)))*ρ) + (k^2) / (Differential(ρ)(V(k, ρ)) + k^2))]
Boundary Conditions: Equation[V(-400000.0, ρ) ~ 2ρ, V(k, 0) ~ 0, Differential(ρ)(V(k, 0)) ~ 2]
Domain: Symbolics.VarDomainPairing[Symbolics.VarDomainPairing(k, -400000.0 .. 300000.0), Symbolics.VarDomainPairing(ρ, 0.0 .. 20000.0)]
Dependent Variables: Num[V(k, ρ)]
Independent Variables: Num[k, ρ]
Parameters: SciMLBase.NullParameters()
Default Parameter ValuesDict{Any, Any}()
julia> NN = 32
32
julia> disc = MOLFiniteDifference([ρ=>NN],k,approx_order=order)
MOLFiniteDifference{MethodOfLines.CenterAlignedGrid, MethodOfLines.ScalarizedDiscretization}(Dict{Num, Int64}(ρ => 32), k, 4, UpwindScheme(1), MethodOfLines.CenterAlignedGrid(), true, false, MethodOfLines.ScalarizedDiscretization(), true, Any[], Base.Pairs{Symbol, Union{}, Tuple{}, @NamedTuple{}}())
julia> prob = discretize(pdesys,disc)
ODEProblem with uType Vector{Float64} and tType Float64. In-place: true
Non-trivial mass matrix: false
timespan: (-400000.0, 300000.0)
u0: 30-element Vector{Float64}:
2580.6451612903224
3870.967741935484
5161.290322580645
6451.612903225807
7741.935483870968
9032.258064516129
10322.58064516129
11612.90322580645
12903.225806451614
14193.548387096775
15483.870967741936
16774.1935483871
18064.516129032258
19354.83870967742
20645.16129032258
21935.483870967742
23225.8064516129
24516.129032258064
25806.451612903227
27096.774193548386
28387.09677419355
29677.41935483871
30967.74193548387
32258.064516129034
33548.3870967742
34838.709677419356
36129.032258064515
37419.354838709674
38709.67741935484
40000.0
julia> sol = solve(prob,Tsit5(),saveat=1e2)
PDESolution:
Return Code:
Success
Dependent variables:
V(k, ρ): (7001, 32) sized solution
Domain:
k ∈ (-400000.0, 300000.0) with 7001 points, step size 100.0
ρ ∈ (0.0, 20000.0) with 32 points, step size 645.1612903225806
From system:
Equations:
L"\begin{align}
k \frac{\mathrm{d}}{\mathrm{d}k} V\left( k, \rho \right) - k^{3} \left( \frac{k^{2}}{\frac{\mathrm{d}}{\mathrm{d}\rho} V\left( k, \rho \right) + k^{2} + 2 \frac{\mathrm{d}}{\mathrm{d}\rho} \frac{\mathrm{d}}{\mathrm{d}\rho} V\left( k, \rho \right) \rho} + \frac{k^{2}}{\frac{\mathrm{d}}{\mathrm{d}\rho} V\left( k, \rho \right) + k^{2}} \right) &= 0
\end{align}
"
Boundary/Initial Conditions:
L"\begin{align}
V\left( -4 \cdot 10^{5}, \rho \right) &= 2 \rho \\
V\left( k, 0 \right) &= 0 \\
\frac{\mathrm{d}}{\mathrm{d}\rho} V\left( k, 0 \right) &= 2
\end{align}
"
julia> plot(sol.disc_data.discretespace.grid[ρ],sol.u[V(k, ρ)][1,:])
julia> plot!(sol.disc_data.discretespace.grid[ρ],sol.u[V(k, ρ)][2,:])
julia> plot!(sol.disc_data.discretespace.grid[ρ],sol.u[V(k, ρ)][3,:])
julia> plot!(sol.disc_data.discretespace.grid[ρ],sol.u[V(k, ρ)][4,:])
one thing that may be odd is that I can’t make V(0, ρ) ~ something
plain code
using ModelingToolkit,MethodOfLines,DomainSets
@parameters k,ρ
@variables V(..)
Dρ = Differential(ρ)
D2ρ = Dρ^2
Dk = Differential(k)
kmin = -4e5
kmax= 3e5
rhomin = 0
rhomax = 2e4
T= 1
N = 2
eq = [k*Dk(V(k,ρ)) ~ T*k^3*(k^2/(k^2+Dρ(V(k,ρ)) + 2*ρ*D2ρ(V(k,ρ))) + (N-1)*k^2/(k^2 + Dρ(V(k,ρ))))]
bcs = [V(kmin,ρ)~λ*ρ, V(k,0)~0,Dρ(V(k,0)) ~ λ ]
@named pdesys = PDESystem(eq,bcs,domains,[k,ρ],[V(k,ρ)])
NN = 32
disc = MOLFiniteDifference([ρ=>NN],k,approx_order=4)
prob = discretize(pdesys,disc)
sol = solve(prob,Tsit5(),saveat=1e2)
plot(sol.disc_data.discretespace.grid[ρ],sol.u[V(k, ρ)][1,:])
plot!(sol.disc_data.discretespace.grid[ρ],sol.u[V(k, ρ)][2,:])
)
.
