Hello, I´m new user of Julia and I´m trying to solve a system of PDES non lineal with method of lines. To check it I´m using a static solution not dependent on time. The idea is to recover to t_i+1 the same solutions. I am going to enumerate the different difficulties and questions that arose and that I did not find answers in the tutorials and pages about the package.
The system is of second order passed to first order. The unknowns are Φ, Λ, Ψ and γ. With
∂ₜ Φ = Φₜ , ∂ₜΨ = Ψₜ , ∂ₜΛ = Λₜ
The code is
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, Plots, DomainSets
#static analytical solutions
Λ_est = x -> 2/(2-(4/5)*x^2*(5-3*x^2))
Φ_est = x -> -(4/5 - 1)*(-2*4/5 + 2)^(1/4)*exp.(-5/2*sqrt(4/5)*atan.((6*x^2 - 5)*sqrt.(4/5)*sqrt.(-25*4/5 + 24)/(25*4/5 - 24))/sqrt.(-25*4/5 + 24) + 5/2*sqrt(4/5)*atan.(sqrt.(4/5)*sqrt.(-25*4/5 + 24)/(25*4/5 - 24))/sqrt.(-25*4/5 + 24))/(3*4/5*x^4 - 5*4/5*x^2 + 2)^(1/4)
Ψ_est = x -> 1
γ_est = x -> 1.0/x * (1.0 - (2/(2-(4/5)*x^2*(5-3*x^2)))^2)
@parameters t x
@variables Φ(..) Φₜ(..) Λ(..) Λₜ(..) Ψ(..) Ψₜ(..) γ(..)
Dt = Differential(t)
Dx = Differential(x)
x_min = 0.001 #is not 0 because the system is irregular at the origin
t_min = 0.0
x_max = 1.0
t_max = 0.1
# q is know
q = x -> (2/5)*(5.0-3.0*x^2)
eq = [Dx(Φ(t,x)) ~ (q(x)*x*Λ(t,x)*Φ(t,x))/(Ψ(t,x)^2),
Dx(Φₜ(t,x)) ~ (q(x)*x)/(Ψ(t,x)^3) * (Φ(t,x)*Ψ(t,x)*Λₜ(t,x) + Φₜ(t,x)*Ψ(t,x)*Λ(t,x) - 2 * Φ(t,x)*Ψₜ(t,x)*Λ(t,x)),
Dt(Ψ(t,x)) ~ Ψₜ(t,x),
Dt(Ψₜ(t,x)) ~ (Φ(t,x)^2/Λ(t,x)^2) * (Dx(Ψ(t,x))/x + (1/2 * γ(t,x) * Ψ(t,x))/x) - ((Ψₜ(t,x)^2)/(2*Ψ(t,x))) + Ψₜ(t,x)*Φₜ(t,x)/Φ(t,x) + Φ(t,x)^2 * ((Dx(Ψ(t,x)))^2/(2*Λ(t,x)*Ψ(t,x)) + q(x)*x*Dx(Ψ(t,x))/(Λ(t,x)*Ψ(t,x)^2) + q(x)/(Λ(t,x)*Ψ(t,x)) + q(x)^2*x^2 /(2*Ψ(t,x)^3)),
Dt(Λ(t,x)) ~ Λₜ(t,x),
Dt(Λₜ(t,x)) ~ -(Φ(t,x)^2/Λ(t,x))*(2*Dx(Ψ(t,x))/(x*Ψ(t,x)) + γ(t,x)/x) + Λ(t,x) * Ψₜ(t,x)/(Ψ(t,x)^2) - Λₜ(t,x) * Φₜ(t,x)/Φ(t,x) - (Φ(t,x)^2/Ψ(t,x)^2)* ( Dx(Ψ(t,x))^2/Λ(t,x) + 2*q(x)*x*Dx(Ψ(t,x))/Ψ(t,x) + 2*q(x) + q(x)^2 * x^2 * Λ(t,x)/(Ψ(t,x)^2)),
Dt(γ(t,x)) ~ (2*Λ(t,x)/Ψ(t,x)^3) *(Ψ(t,x)^3*Λₜ(t,x)+ q(x)*x*Λ(t,x)^2*Ψₜ(t,x) - Dx(Ψₜ(t,x)))
]
bcs=[
Ψ(t_min,x) ~ 1.0,
Ψₜ(t_min,x) ~ 0.0,
Λ(t_min,x) ~ 2/(2-(4/5)*x^2*(5-3*x^2)),
Λₜ(t_min,x) ~ 0.0,
γ(t_min,x) ~ (1.0/x) * (1.0 - (2/(2-(4/5)*x^2*(5-3*x^2)))^2),
Φ(t,x_min) ~ Φ_est(x_min),
Φₜ(t,x_min) ~ 0.0,
Φₜ(t_min,x) ~ 0.0,
Φ(t_min,x) ~ -(4/5 - 1)*(-2*4/5 + 2)^(1/4)*exp.(-5/2*sqrt.(4/5)*atan.((6*x^2 - 5)*sqrt.(4/5)*sqrt.(-25*4/5 + 24)/(25*4/5 - 24))/sqrt.(-25*4/5 + 24) + 5/2*sqrt(4/5)*atan.(sqrt.(4/5)*sqrt.(-25*4/5 + 24)/(25*4/5 - 24))/sqrt.(-25*4/5 + 24))/(3*4/5*x^4 - 5*4/5*x^2 + 2)^(1/4)
]
domains = [t ∈Interval(t_min, t_max),
x ∈Interval(x_min, x_max)]
@named pdesys = PDESystem(eq, bcs, domains,[t, x],[Φ(t,x), Φₜ(t,x), Ψ(t,x), Ψₜ(t,x), Λ(t,x), Λₜ(t,x),γ(t,x)])
N = 32
dx = (x_max -x_min)/N # intervalo de discretización
discx = chebyspace(N, domains[2])
discretization = MOLFiniteDifference([x => N],t)
@time prob = discretize(pdesys, discretization)
alg = Rosenbrock23()
@time sol = solve(prob,alg)
discrete_x = sol[x]
discrete_t = sol[t]
solΦ = sol[Φ(t, x)]
solΦₜ = sol[Φₜ(t, x)]
solΛ = sol[Λ(t, x)]
solΛₜ = sol[Λₜ(t, x)]
solΨ = sol[Ψ(t, x)]
solΨₜ = sol[Ψₜ(t, x)]
Φ_real = [Φ_est(s) for s in discrete_x]
errΦ = scatter()
for i in 1:length(discrete_t)
scatter!(discrete_x, solΦ[i, :].-Φ_real,title= "Error en Φ", label="t=$(discrete_t[i])")
end
display(errΦ)
The code runs, one of the problems I have is that if I want to overcome the barrier of N greater than 100 I have memory allocation problems and it crashes. Do you know why this happens?. Is it because of the complexity of the system? Or is there something that as a newbie I don’t do (I must clarify that I do not have a great computer).
On the other hand, if I try to use a smaller grid but with a chebyshev-type discretization Non-Uniform Rectilinear Grids · MethodOfLines.jl , when I run the pdes discretization, I get the following error:
BoundsError: attempt to access 32-element Vector{StaticArraysCore.SVector{2, Float64}} at index [33]
Lastly, as can be seen, I do not have an evolution equation for Φ, so errors accumulate in its calculation and it does not give a static solution. My question is why when I give only the bc Φ(x=0.001,t) = Φ_est(0.001) and Φ_t(x=0.001,t) = 0 it does not solve for Phi and Phi_t ? . But if I give it those conditions plus Φ(x,t=0) = Φ_est and Φ(x,t=0) = 0.0 if it gives me valid results
I would appreciate if at least one doubt can be resolved in advance
Sorry if my writing is bad, english is not my native language, thank you very much