ExtraVariableSystemException while working with DiffEqOperators.jl

Specific Domains Modelling & Simulations
1

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[2](τ)
f[3](τ)
#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.

Thanks in advance for your help.

1 Like
2

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 :sweat_smile: . It will be really cool in 2023 though.

2 Likes
3

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

1 Like