# ExtraVariableSystemException while working with DiffEqOperators.jl

I’ve run into two issues while using DiffEqOperators.jl to solve the Feynman-Kac Formula through the method of lines.

Here is the setup to the problem:

``````#setup symbols
@parameters τ y
@variables f(..)
Dτ = Differential(τ)
Dy = Differential(y)
Dyy = Differential(y)^2

#setup parameters

ρ = 0.025

#parameters
σ = 5.0
θ = 3.0
μ = 50.0

λ = 0.5

#bounds
T = 20
y_min = 1e4
y_max = 1e4 +20

#domains
dom = [ τ ∈ Interval(0.0,T),
y ∈ Interval(y_min,y_max)]

#define the feynman-kac equation
equation = Dτ(f(τ,y)) + θ*(μ-y)*Dy(f(τ,y)) + (σ^2)/2 * Dyy(f(τ,y)) ~ - ρ*f(τ,y)

boundary_conditions = [
f(T,y) ~ 0.6
,f(τ,y_min) ~ 0 #no probability at minimum
,f(τ,y_max) ~ 0 #No probability at the high point
]

fk = PDESystem(equation,boundary_conditions,dom,[τ,y],[f(τ,y)],name=:feynman_kac);

Δy = 1 #just to get the error to show up easily
discrete = MOLFiniteDifference([y=>Δy]
,τ
;centered_order=2
)
``````

when I attempt to discretize the problem using

``````discretize(fk,discrete)
``````

I get

``````ExtraVariablesSystemException: The system is unbalanced. There are 21 highest order derivative variables and 19 equations.
More variables than equations, here are the potential extra variable(s):
f(τ)
f(τ)
#etc
``````

Depending on the specific choice of boundaries and \Delta y I’ll get either an off by one or off by two error on the derivative variables vs equations.

I’m understanding that there is a difference in the dimensions of the estimated derivatives and the equation approximation. Is that correct?

I suspect it is some oversight on my part as I’m quite new to solving PDEs and julia in general.

1 Like

It’s not you, the discretizer is just a bit dumb right now. There’s a big warning sign on the top of the documentation that it’s still heavily being worked on, and that will be removed when it’s stable. Until then, there’s a lot of edge cases.

Here, the problem is that it currently doesn’t reorder variables. So you need to define the equation as:

``````equation = Dτ(f(τ,y)) ~ - ρ*f(τ,y) - θ*(μ-y)*Dy(f(τ,y)) - (σ^2)/2 * Dyy(f(τ,y))
``````

so it picks up that it’s an ODE. But then you have to flip the signs because it doesn’t recognize the terminal condition at `T`, and wants an initial condition. So this works:

``````using DiffEqOperators, ModelingToolkit, DomainSets

#setup symbols
@parameters τ y
@variables f(..)
Dτ = Differential(τ)
Dy = Differential(y)
Dyy = Differential(y)^2

#setup parameters

ρ = 0.025

#parameters
σ = 5.0
θ = 3.0
μ = 50.0

λ = 0.5

#bounds
T = 20
y_min = 1e4
y_max = 1e4 +20

#domains
dom = [ τ ∈ Interval(0.0,T),
y ∈ Interval(y_min,y_max)]

#define the feynman-kac equation
equation = Dτ(f(τ,y)) ~ ρ*f(τ,y) + θ*(μ-y)*Dy(f(τ,y)) + (σ^2)/2 * Dyy(f(τ,y))

boundary_conditions = [
f(0.0,y) ~ 0.6
,f(τ,y_min) ~ 0 #no probability at minimum
,f(τ,y_max) ~ 0 #No probability at the high point
]

fk = PDESystem(equation,boundary_conditions,dom,[τ,y],[f(τ,y)],name=:feynman_kac);

Δy = 1 #just to get the error to show up easily
discrete = MOLFiniteDifference([y=>Δy]
,τ
;centered_order=2
)
discretize(fk,discrete)
``````

and that explains the reason why the documentation said it’s still being developed so use at your own risk . It will be really cool in 2023 though.

Thanks for the follow up, I’m giving it a try.