# Trying to use data as an initial condition for a 1D diffusion problem with DiffEqOperators.jl

Using code from Example 1 in the library repo README to get started.

I’m getting `ERROR: MethodError: no method matching zero(::Type{Vector{Float64}})` and not clear how to decode the Stacktrace …

``````using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets, DifferentialEquations, Interpolations
###
function uniform_linear_sample(data_in, n_samples)
# Resample a 1D dataset using a given number of sample points using linear interpolation.

nodes_in = range(1, size(data_in)[1], length=length(data_in))
nodes_eval = range(1, Int(nodes_in[end]), length = Int(n_samples))

int = interpolate((nodes_in,), data_in, Gridded(Linear()))

data_out = zeros(length(nodes_eval))
idx = trues(length(nodes_eval))

data_out[idx] = int(nodes_eval[idx])

return data_out
end

n_nodes = 1000

x_ic = range(0,Float64(2pi),length=10000)
data_vector = cos.(x_ic) .+ 1

ic() = uniform_linear_sample(data_vector,n_nodes)

@register ic()
###
# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
bcs = [u(0,x) ~ ic(),
u(t,0) ~ exp(-t),
u(t,Float64(pi)) ~ -exp(-t)]

# Space and time domains
domains = [t ∈ Interval(0.0,1.0),
x ∈ Interval(0.0,Float64(2pi))]

# PDE system
@named pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x=>dx],t;centered_order=order)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)

# Solve ODE problem
sol = solve(prob,saveat=0.1)
``````
``````ERROR: LoadError: MethodError: no method matching zero(::Type{Vector{Float64}})
Closest candidates are:
zero(::Union{Type{P}, P}) where P<:Dates.Period at C:\buildbot\worker\package_win64\build\usr\share\julia\stdlib\v1.6\Dates\src\periods.jl:53
zero(::ExtendedJumpArray) at C:\Users\thompsmj.X-MJIS3050PC09\.julia\packages\DiffEqJump\x05Qi\src\extended_jump_array.jl:26
zero(::Unitful.AbstractQuantity{T, D, U} where {T, D, U<:(Unitful.Units{N, D, A} where {N, D, A<:Unitful.Affine})}) at C:\Users\thompsmj.X-MJIS3050PC09\.julia\packages\Unitful\yI5QN\src\quantities.jl:374
...
Stacktrace:
[1] zero(x::Vector{Vector{Float64}})
@ Base .\abstractarray.jl:1085
``````

make it `ic(x)` and use a callable.

Hi Chris, thanks for trying to help, but I’m still unclear. I can interpolate the data fine, but I don’t understand what to do to “use a callable,” though I’ve been going through the docs to try to get the hang of it.

I tried using exactly what I have posted above, and swapping in this line as you suggest:

``````bcs = [u(0,x) ~ ic(x),
u(t,0) ~ exp(-t),
u(t,Float64(pi)) ~ -exp(-t)]
``````

And I get

``````ERROR: LoadError: MethodError: no method matching ic(::Num)
``````

on that line.

Above that, I then tried

``````ic(x) = uniform_linear_sample(data_vector,n_nodes)
``````

which gets through to the end and gives this error at the `solve` line

``````sol = solve(prob,saveat=0.1)

ERROR: LoadError: MethodError: no method matching zero(::Type{Vector{Float64}})
``````

Edit: fixed link to wrong place.

Why not just use the DataInterpolations.jl linear interpolation?

Didn’t know about it! Looks like a useful library.

Trying it out though I’m getting `LoadError: UndefVarError: LinearInterpolation not defined` with the current v3.6.1.

Show what you ran.

OK, the issue with `LinearInterpolation not defined` is resolved if I don’t load `Interpolations` and `DataInterpolations` together.

``````using DataInterpolations, Interpolations
u = rand(5)
t = 0:4
interp = LinearInterpolation(u,t)

ERROR: LoadError: UndefVarError: LinearInterpolation not defined
``````
``````using DataInterpolations
u = rand(5)
t = 0:4
interp = LinearInterpolation(u,t)

10-element LinearInterpolation{Vector{Float64}, UnitRange{Int64}, true, Float64}:
``````

This is where I’m unfamiliar with symbolic programming, so the final line gives `LoadError: MethodError: no method matching munge_data(::Vector{Float64}, ::Float64)`

running:

``````using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets, DifferentialEquations, DataInterpolations

x_dense = range(0,Float64(2pi),length=10000)
data = cos.(x_dense) .+ 1.0 .+ 0.05*randn(10000)

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
ic(x) = LinearInterpolation(data, x)
@register ic(x)
eq  = Dt(u(t,x)) ~ Dxx(u(t,x))
bcs = [u(0,x) ~ ic(x),
u(t,0) ~ exp(-t),
u(t,Float64(pi)) ~ -exp(-t)]

# Space and time domains
domains = [t ∈ Interval(0.0,1.0),
x ∈ Interval(0.0,Float64(2pi))]

# PDE system
@named pdesys = PDESystem(eq,bcs,domains,[t,x],[u(t,x)])

# Method of lines discretization
dx = 0.1
order = 2
discretization = MOLFiniteDifference([x=>dx],t;centered_order=order)

# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)
``````

Also tried `ic(x) = LinearInterpolation(data, 0:Float64(2pi))` with all else the same, and running `prob = discretize(pdesys,discretization)` succeeds, giving:

``````ODEProblem with uType Vector{LinearInterpolation{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, true, Float64}} and tType Float64. In-place: true
``````

but `sol = solve(prob,Tsit5(),saveat=0.1)` gives:

``````ERROR: LoadError: MethodError: Cannot `convert` an object of type Vector{Float64} to an object of type LinearInterpolation{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, true, Float64}
``````

No, LinearInterpolation is a callable. Look at the documentation.

``````using DataInterpolations
u = rand(5)
t = 0:4
interp = LinearInterpolation(u,t)
interp(t)
``````

Because of that, you do:

``````ic = LinearInterpolation(data, x_dense)
``````

now `ic(x)` is the interpolation function over the data, which means

``````@register ic(x)
bcs = [u(0,x) ~ ic(x),
u(t,0) ~ exp(-t),
u(t,Float64(pi)) ~ -exp(-t)]
``````

will work.

``````using OrdinaryDiffEq, ModelingToolkit, DiffEqOperators, DomainSets, DifferentialEquations, DataInterpolations

x_dense = range(0,Float64(2pi),length=10000)
data = cos.(x_dense) .+ 1.0

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)
Dt = Differential(t)
Dxx = Differential(x)^2

# 1D PDE and boundary conditions
ic = LinearInterpolation(data, x_dense)
@register ic(x)
``````
``````ERROR: LoadError: cannot define function ic; it already has a value
``````

`ic` iteself is:

``````julia> ic
20000-element LinearInterpolation{Vector{Float64}, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}}, true, Float64}:
``````

`ic(x)` is:

``````julia> ic(x)
ERROR: MethodError: isless(::Num, ::Float64) is ambiguous. Candidates:
isless(x::Real, y::AbstractFloat) in Base at operators.jl:168
isless(a::Num, b::Real) in Symbolics at C:\Users\thompsmj.X-MJIS3050PC09\.julia\packages\Symbolics\Cmx10\src\num.jl:130
Possible fix, define
isless(::Num, ::AbstractFloat)
``````
``````import Symbolics
(interp::DataInterpolations.AbstractInterpolation)(t::Symbolics.Num) = Symbolics.SymbolicUtils.term(interp,t)
``````

and the release with that will make it automatic.

I’ve been running into similar issues using DataInterpolations’s `CubicSplines` in ModelingToolkit. The pull request provided by @ChrisRackauckas looks like it may fix my problems, but I can’t update Symbolics to the latest version. There’s some sort of incompatibility with the packages in ModelingToolkit. Any thoughts?

``````pkg> add Symbolics@4.0
Resolving package versions...
ERROR: Unsatisfiable requirements detected for package ModelingToolkit [961ee093]:
ModelingToolkit [961ee093] log:
├─possible versions are: 0.0.1-6.7.1 or uninstalled
├─restricted to versions * by QED [03dc1bc1], leaving only versions 0.0.1-6.7.1
│ └─QED [03dc1bc1] log:
│   ├─possible versions are: 0.1.0 or uninstalled
│   └─QED [03dc1bc1] is fixed to version 0.1.0
├─restricted by compatibility requirements with Symbolics [0c5d862f] to versions: 0.0.1-5.9.1 or uninstalled, leaving only versions: 0.0.1-5.9.1
│ └─Symbolics [0c5d862f] log:
│   ├─possible versions are: 0.1.0-4.0.0 or uninstalled
│   ├─restricted to versions * by QED [03dc1bc1], leaving only versions 0.1.0-4.0.0
│   │ └─QED [03dc1bc1] log: see above
│   └─restricted to versions 4.0 by an explicit requirement, leaving only versions 4.0.0
├─restricted by compatibility requirements with RuntimeGeneratedFunctions [7e49a35a] to versions: [0.0.1-3.21.0, 4.5.0-6.7.1] or uninstalled, leaving only versions: [0.0.1-3.21.0, 4.5.0-5.9.1]
│ └─RuntimeGeneratedFunctions [7e49a35a] log:
│   ├─possible versions are: 0.1.0-0.5.3 or uninstalled
│   └─restricted by compatibility requirements with Symbolics [0c5d862f] to versions: 0.5.0-0.5.3
│     └─Symbolics [0c5d862f] log: see above
├─restricted by compatibility requirements with Latexify [23fbe1c1] to versions: [0.0.1-0.8.0, 1.2.7-6.7.1] or uninstalled, leaving only versions: [0.0.1-0.8.0, 1.2.7-3.21.0, 4.5.0-5.9.1]
│ └─Latexify [23fbe1c1] log:
│   ├─possible versions are: 0.5.0-0.15.9 or uninstalled
│   ├─restricted by compatibility requirements with Symbolics [0c5d862f] to versions: 0.11.0-0.15.9
│   │ └─Symbolics [0c5d862f] log: see above
│   └─restricted by compatibility requirements with Requires [ae029012] to versions: 0.13.0-0.15.9 or uninstalled, leaving only versions: 0.13.0-0.15.9
│     └─Requires [ae029012] log:
│       ├─possible versions are: 0.5.0-1.1.3 or uninstalled
│       └─restricted by compatibility requirements with Symbolics [0c5d862f] to versions: 1.1.0-1.1.3
│         └─Symbolics [0c5d862f] log: see above
├─restricted by compatibility requirements with StaticArrays [90137ffa] to versions: 4.0.7-6.7.1 or uninstalled, leaving only versions: 4.5.0-5.9.1
│ └─StaticArrays [90137ffa] log:
│   ├─possible versions are: 0.8.0-1.2.13 or uninstalled
│   ├─restricted to versions * by Equilibrium [187daef3], leaving only versions 0.8.0-1.2.13
│   │ └─Equilibrium [187daef3] log:
│   │   ├─possible versions are: 0.1.0 or uninstalled
│   │   └─Equilibrium [187daef3] is fixed to version 0.1.0
│   ├─restricted by compatibility requirements with ForwardDiff [f6369f11] to versions: 0.8.3-1.2.13
│   │ └─ForwardDiff [f6369f11] log:
│   │   ├─possible versions are: 0.9.0-0.10.21 or uninstalled
│   │   ├─restricted to versions * by Equilibrium [187daef3], leaving only versions 0.9.0-0.10.21
│   │   │ └─Equilibrium [187daef3] log: see above
│   │   ├─restricted by compatibility requirements with DiffRules [b552c78f] to versions: 0.10.6-0.10.21 or uninstalled, leaving only versions: 0.10.6-0.10.21
│   │   │ └─DiffRules [b552c78f] log:
│   │   │   ├─possible versions are: 0.0.8-1.3.1 or uninstalled
│   │   │   ├─restricted by compatibility requirements with ForwardDiff [f6369f11] to versions: [0.0.8-1.0.2, 1.2.1-1.3.1]
│   │   │   │ └─ForwardDiff [f6369f11] log: see above
│   │   │   └─restricted by compatibility requirements with Symbolics [0c5d862f] to versions: 0.1.0-1.3.1, leaving only versions: [0.1.0-1.0.2, 1.2.1-1.3.1]
│   │   │     └─Symbolics [0c5d862f] log: see above
│   │   └─restricted by compatibility requirements with StaticArrays [90137ffa] to versions: 0.10.13-0.10.21 or uninstalled, leaving only versions: 0.10.13-0.10.21
│   │     └─StaticArrays [90137ffa] log: see above
│   └─restricted by compatibility requirements with Symbolics [0c5d862f] to versions: 1.1.0-1.2.13
│     └─Symbolics [0c5d862f] log: see above
└─restricted by compatibility requirements with SymbolicUtils [d1185830] to versions: 0.0.1-3.1.1 or uninstalled — no versions left
└─SymbolicUtils [d1185830] log:
├─possible versions are: 0.1.0-0.18.0 or uninstalled
└─restricted by compatibility requirements with Symbolics [0c5d862f] to versions: 0.18.0
└─Symbolics [0c5d862f] log: see above
``````

Symbolics.jl v4 was just released, and this requires that update. So some downstream updates will need to be done before this can really be used everywhere. Working on it

1 Like

Sounds good, thanks!