I would like to solve a PDE with MethodOfLines.jl and ModellingToolkit.jl similar to the example docs.sciml.ai/MethodOfLines/stable/tutorials/PIDE/. My modification is to change the integral term. Instead of integrating u(t,x) wrt x I would like to integrate (y-x)u(t,y) wrt y. Thus, I have to evaluate the function u which is originally defined as a function of x, with a paramter y in the integral. u still also needs to be able to be evaluated with parameter x since in the other parts of the PDE I still have to evaluate u(t,x). How can I do that? Because when I just define y as a paramter and then write u(t,y) I get the error message:
ERROR: ArgumentError: invalid index: nothing of type Nothing
Stacktrace:
[1] to_index(i::Nothing)
@ Base ./indices.jl:300
[2] to_index(A::Vector{Symbolics.VarDomainPairing}, i::Nothing)
@ Base ./indices.jl:277
[3] _to_indices1(A::Vector{Symbolics.VarDomainPairing}, inds::Tuple{Base.OneTo{Int64}}, I1::Nothing)
@ Base ./indices.jl:359
[4] to_indices
@ Base ./indices.jl:354 [inlined]
[5] to_indices
@ Base ./indices.jl:345 [inlined]
[6] getindex
@ Base ./abstractarray.jl:1288 [inlined]
[7] (::PDEBase.var"#54#62"{Vector{Symbolics.VarDomainPairing}})(x::SymbolicUtils.BasicSymbolic{Real})
@ PDEBase ~/.julia/packages/PDEBase/VVGPP/src/variable_map.jl:34
[8] iterate
@ Base ./generator.jl:47 [inlined]
[9] collect_to!(dest::Vector{Pair{…}}, itr::Base.Generator{Vector{…}, PDEBase.var"#54#62"{…}}, offs::Int64, st::Int64)
@ Base ./array.jl:892
[10] collect_to_with_first!(dest::Vector{…}, v1::Pair{…}, itr::Base.Generator{…}, st::Int64)
@ Base ./array.jl:870
[11] _collect(c::Vector{…}, itr::Base.Generator{…}, ::Base.EltypeUnknown, isz::Base.HasShape{…})
@ Base ./array.jl:864
[12] collect_similar(cont::Vector{Any}, itr::Base.Generator{Vector{Any}, PDEBase.var"#54#62"{Vector{…}}})
@ Base ./array.jl:763
[13] map(f::Function, A::Vector{Any})
@ Base ./abstractarray.jl:3282
[14] PDEBase.VariableMap(pdesys::PDESystem, disc::MOLFiniteDifference{…}; replaced_vars::Dict{…})
@ PDEBase ~/.julia/packages/PDEBase/VVGPP/src/variable_map.jl:33
[15] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase ~/.julia/packages/PDEBase/VVGPP/src/symbolic_discretize.jl:17
[16] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
@ PDEBase ~/.julia/packages/PDEBase/VVGPP/src/discretization_state.jl:57
[17] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
@ PDEBase ~/.julia/packages/PDEBase/VVGPP/src/discretization_state.jl:54
[18] top-level scope
@ ~/Documents/implement/multid_op_form/example_1.jl:27
Some type information was truncated. Use show(err) to see complete types.
Here is the code I used:
using MethodOfLines, ModelingToolkit, OrdinaryDiffEq, DomainSets, Plots
@parameters t, x, y
@variables u(..) cumuSum(..) integrand(..)
Dt = Differential(t)
Dx = Differential(x)
xmin = 0.0
xmax = 2.0pi
# Integral limits are defined with DomainSets.ClosedInterval
Iy = Integral(y in DomainSets.ClosedInterval(xmin, xmax)) # basically cumulative sum from 0 to x
eq = [
integrand(t,x, y) ~ (y-x)*u(t, y), # Define the integrand as a separate function
cumuSum(t, x) ~ Iy(integrand(t,x,y)), # Note wrapping the argument to the derivative with an auxiliary variable
Dt(u(t, x)) + 2 * u(t, x) + 5 * Dx(cumuSum(t, x)) ~ 1
]
bcs = [u(0.0, x) ~ cos(x), Dx(u(t, xmin)) ~ 0.0, Dx(u(t, xmax)) ~ 0]
domains = [t ∈ Interval(0.0, 2.0), x ∈ Interval(xmin, xmax)]
@named pde_system = PDESystem(eq, bcs, domains, [t, x], [u(t, x), cumuSum(t, x)])
order = 2
discretization = MOLFiniteDifference([x => 30], t)
prob = MethodOfLines.discretize(pde_system, discretization)
sol = solve(prob, QNDF(), saveat = 0.1);
solu = sol[u(t, x)]
plot(sol[x], transpose(solu))
