Comments (5)
Transferring to MethodOfLines as this is more of a MethodOfLines issue. If it generates the right DAE then MTK should be fine.
from methodoflines.jl.
I have the same Problem with another system.
Discretize complains that there are no valid boundary conditions
@chris
Even PDAEs can be very tricky, it would be good if we had a simple example that works
from methodoflines.jl.
I came across a simple PDAE problem using Merton's portfolio problem (without consumption and discounting). Since the boundary can be solved explicitly, discretize
can easily chop the PDAE into a DAE.
However, all the solvers that I've tried (Rosenbrock23
, QNDF
, Rodas4P
, Rodas5P
) fail with MaxIters
during the first time step. Perhaps I am doing something wrong, but I just can't spot the mistake.
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Plots
@parameters x t
@variables J(..) m(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
r::Float64 = 0.03 # bond rate
μ::Float64 = r + 0.07 # expected market return
σ::Float64 = 0.2 # market volatility
γ::Float64 = 2 # risk aversion
x_min::Float64 = 5
x_max::Float64 = 100
t_min::Float64 = 0.0
t_max::Float64 = 40
# Notice that time is reversed
m_analytical(x, t) = (μ - r) / (γ * (σ^2))
J_analytical(x, t) = (exp(t .* (1 .- γ) .* (r + (1 / (2 * γ)) * (((μ - r)^2) / ((σ^2)))))) .* (x .^ (1 .- γ)) ./ (1 .- γ)
JT(x) = (x^(1 - γ)) / (1 - γ)
eq = [
# HJB
0 ~ -Dt(J(x, t)) + Dx(J(x, t)) * (x * (r + m(x, t) * (μ - r))) + (1 / 2) * Dxx(J(x, t)) * (x^2) * ((m(x, t))^2) * (σ^2),
# FOC
0 ~ Dx(J(x, t)) * x * (μ - r) + Dxx(J(x, t)) * (x^2) * (m(x, t)) * (σ^2),
]
domains = [x ∈ Interval(x_min, x_max),
t ∈ Interval(t_min, t_max)]
# Boundary conditions
bcs = [
J(x, t_min) ~ JT(x),
m(x, t_min) ~ (μ - r) / (γ * (σ^2))
]
@named pdesys = PDESystem(eq, bcs, domains, [x, t], [J(x, t), m(x, t)])
#
# Method of lines discretization
#
Nx::Int64 = 32
order::Int64 = 2 # This may be increased to improve accuracy of some schemes
# Use a Float to directtly specify stepsizes dx.
discretization = MOLFiniteDifference([x => Nx], t, approx_order=order)
# Convert the PDAE problem into an DAE problem
println("Discretization:")
@time prob = discretize(pdesys, discretization)
#
# Solving the problem
#
println("Solve:")
@time sol = solve(prob, Rodas4P(), saveat=0.1)
from methodoflines.jl.
Is the system structurally simplified? Is there a Jacobian singularity at the start?
from methodoflines.jl.
First of all, I want to thank you @ChrisRackauckas for your time, and warn you that I'm a complete noob when it comes to solving PDEs. So feel free to close the issue if you believe that it becomes sterile.
A singularity seems quite likely indeed, since
PDAE formulation code
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Sundials, ODEInterfaceDiffEq, LSODA, ForwardDiff, Plots
@parameters x t
@variables J(..) m(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
bond_rate::Float64 = 0.03
stock_expected_return::Float64 = bond_rate + 0.07
stock_volatility::Float64 = 0.2
risk_aversion::Float64 = 2
initial_wealth::Float64 = 10
external_portfolio_contribution::Float64 = 8
x_min::Float64 = 5
x_max::Float64 = 100
t_min::Float64 = 0.0
t_max::Float64 = 40
JT(x) = (x^(1 - risk_aversion)) / (1 - risk_aversion)
analytical_A = ((-(1 - risk_aversion) .* (bond_rate + (((stock_expected_return - bond_rate) ./ stock_volatility) .^ 2) ./ (2 .* risk_aversion))) ./ risk_aversion)
analytical_g(t) = exp.(.-analytical_A .* t)
analytical_J(x, t) = ((analytical_g(t)) .^ (risk_aversion)) .* (((x .+ external_portfolio_contribution .* (1 .- exp.(-t .* bond_rate)) ./ (bond_rate)) .^ (1 .- risk_aversion)) ./ (1 .- risk_aversion))
analytical_m(x, t) = ((stock_expected_return .- bond_rate) ./ (stock_volatility .^ 2)) .* ((x .+ external_portfolio_contribution .* (1 .- exp.(.-t .* bond_rate)) / (bond_rate)) ./ (x .* risk_aversion))
eq = [
# HJB
0 ~ -Dt(J(x, t)) + Dx(J(x, t)) * (x * (bond_rate + m(x, t) * (stock_expected_return - bond_rate)) + external_portfolio_contribution) + (1 / 2) * Dxx(J(x, t)) * (x^2) * ((m(x, t))^2) * (stock_volatility^2),
# FOC
0 ~ Dx(J(x, t)) * x * (stock_expected_return - bond_rate) + Dxx(J(x, t)) * (x^2) * (m(x, t)) * (stock_volatility^2),
]
domains = [x ∈ Interval(x_min, x_max),
t ∈ Interval(t_min, t_max)]
# Boundary conditions
bcs = [
J(x, t_min) ~ JT(x),
m(x, t_min) ~ (stock_expected_return - bond_rate) / (risk_aversion * (stock_volatility^2))
]
@named pdesys = PDESystem(eq, bcs, domains, [x, t], [J(x, t), m(x, t)])
#
# Method of lines discretization
#
Nx::Int64 = 128
order::Int64 = 2 # This may be increased to improve accuracy of some schemes
# Use a Float to directtly specify stepsizes dx.
# discretization = MOLFiniteDifference([x => Nx], t, approx_order=order)
discretization = MOLFiniteDifference([x => Nx], t, approx_order=order, use_ODAE=true)
# Convert the PDE problem into an ODE problem
println("Discretization:")
@time prob = discretize(pdesys, discretization)
#
# Solving the problem
#
println("Solve:")
@time sol = solve(prob, Sundials.IDA(), saveat=0.1)
Another problem with the following PDAE formulation above is that MOLFiniteDifference(..., use_ODAE=true)
produces an ODESystem
instead of a DAE (I assume it should be a DAE because Sundials.IDA()
throws an error.
The system has an equivalent PDE formulation whenever Sundials
and ODEInterfaceDiffEq
solvers can solve 14 years, while native julia solvers only manage to find the solution for the first 5 years.
The problem with this equivalent PDE formulation is that I'm unable to reliably extract
Equivalent PDE formulation code
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, Sundials, ODEInterfaceDiffEq, LSODA, ForwardDiff, FiniteDifferences, Plots
@parameters x t
@variables J(..)
Dt = Differential(t)
Dx = Differential(x)
Dxx = Differential(x)^2
bond_rate::Float64 = 0.03
stock_expected_return::Float64 = bond_rate + 0.07
stock_volatility::Float64 = 0.2
risk_aversion::Float64 = 2
initial_wealth::Float64 = 10
external_portfolio_contribution::Float64 = 8
x_min::Float64 = 5
x_max::Float64 = 100
t_min::Float64 = 0.0
t_max::Float64 = 40
# t_max::Float64 = 10
JT(x) = (x^(1 - risk_aversion)) / (1 - risk_aversion)
analytical_A = ((-(1 - risk_aversion) .* (bond_rate + (((stock_expected_return - bond_rate) ./ stock_volatility) .^ 2) ./ (2 .* risk_aversion))) ./ risk_aversion)
analytical_g(t) = exp.(.-analytical_A .* t)
analytical_J(x, t) = ((analytical_g(t)) .^ (risk_aversion)) .* (((x .+ external_portfolio_contribution .* (1 .- exp.(-t .* bond_rate)) ./ (bond_rate)) .^ (1 .- risk_aversion)) ./ (1 .- risk_aversion))
eq = [
# HJB
0 ~ -Dt(J(x, t)) + Dx(J(x, t)) * (x * bond_rate + external_portfolio_contribution) - (1 / 2) * ((stock_expected_return - bond_rate)^2) / (stock_volatility^2) * ((Dx(J(x, t))^2) / Dxx(J(x, t))),
]
domains = [x ∈ Interval(x_min, x_max),
t ∈ Interval(t_min, t_max)]
# Boundary conditions
bcs = [
J(x, t_min) ~ JT(x)
]
@named pdesys = PDESystem(eq, bcs, domains, [x, t], [J(x, t)])
#
# Method of lines discretization
#
Nx::Int64 = 128
order::Int64 = 2 # This may be increased to improve accuracy of some schemes
# Use a Float to directtly specify stepsizes dx.
# discretization = MOLFiniteDifference([x => Nx], t, approx_order=order)
discretization = MOLFiniteDifference([x => Nx], t, approx_order=order, jac=true) # It seems that jac is not useful at all
# Convert the PDE problem into an ODE problem
println("Discretization:")
# @time prob = discretize(pdesys, discretization)
@time prob = discretize(pdesys, discretization, jac=true) # It seems that jac is not useful at all
#
# Solving the problem
#
println("Solve:")
@time sol = solve(prob, Sundials.CVODE_BDF(), saveat=0.1)
# The solution is actually quite accurate
# It seems that interpolation doesn't include derivatives, since this command fails
sol(20.0, 8.0, deriv=Val{1})
# Trying to obtain the second derivative Jxx using ForwardDiff gives 0, which is wrong
ForwardDiff.derivative(x_val -> ForwardDiff.derivative(x_val -> sol(x_val, 8.0, dv=J(x, t)), x_val), 20.0)
# Using finite differences to approximate the nonconstant component of m(x,t) yields nonsense
forward_fdm(12, 1, max_range=9e-20)(x_val -> sol(x_val, 8.0, dv=J(x, t)), 20.0) / (- 20 * forward_fdm(12, 2, max_range=9e-20)(x_val -> sol(x_val, 8.0, dv=J(x, t)), 20.0))
from methodoflines.jl.
Related Issues (20)
- Interpolate grid values for staggered grid discretization
- Symbolic tracing in staggered grid can return Nans
- Move to using SymbolicUtils chains
- Package installation error HOT 1
- `discretize` errrors when a subset of equations have no time derivatives HOT 5
- Equation + State mismatch dependent on grid spacing
- MethodOfLines down-versions ModelingToolkit HOT 1
- Error with test/pde_systems/MOL_1D_Interface_Nonlinear.jl
- MOL & BoundaryConditions HOT 3
- PDESystem & arguments
- "ExtraVariablesSystemException: The system is unbalanced" Error while solving a PDE! HOT 4
- System of PDEs Initial Failure warning and MethodError HOT 1
- System of PDEs `solve` error
- remake problem tutorial lets user accidentally re-order parameters HOT 3
- Complex boundary conditions HOT 1
- MoL : AssertionError: Boundary condition is not on a boundary of the domain, or is not a valid boundary condition
- MethodOfLines issue with coupled non-linear PDE system
- AssertionError: Boundary condition (...) is not on a boundary of the domain, or is not a valid boundary condition
- ArgumentError: Differential w.r.t. multiple variables Any[0.030,....] are not allowed.
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from methodoflines.jl.